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Effect of Uniformly Distributed Roughness on 
Turbulent Skin-Friction Drag at Supersonic 
Speeds' 


FRANK E. GODDARD, JR.** 
California Institute of Technology 


SUMMARY 


An experimental program was carried out in the 18-in. by 20- 
in. supersonic wind tunnel of the Jet Propulsion Laboratory to 
determine the effect of uniformly distributed sand-grain rough- 
ness on the skin-friction drag of a body of revolution for the case 
of a turbulent boundary layer. The Mach Number range cov- 
ered was 1.98 to 4.54, and the Reynolds Number varied from 
about 3 X 108 to 8 X 108 Some data were also obtained at a 
Mach Number of 0.70. 

At speeds up to a Mach Number of 5 and for roughness sizes 
such that the quadratic resistance law holds, the compressibility 
effect is indirect, and the skin-friction drag is a function of only 
the roughness Reynolds Number vsk/v, exactly as in the incom- 
pressible case. It is shown that the entire compressibility effect 
is a reduction of the fluid density at the surface as the Mach 
Number increases. The critical roughness, below which the 
surface is hydraulically smooth, is &eriticas ~ 10v/v% ; this is equal 
to the thickness of the laminar sublayer for a smooth surface for 
both compressible and incompressible flow. Over the range of 
roughness sizes considered here, there appears to be no wave drag 
associated with the drag due to roughness. The shift in the tur- 
bulent veocity profile A(u/vs ) for a rough surface at supersonic 
speeds is a function of only the roughness Reynolds Number 
v%k/v and quantitatively follows exactly the same law as that 


for the incompressible case. 


SYMBOLS 


afterbody 
drag coefficient of a single sand grain 
average skin-friction-drag coefficient 
local skin-friction-drag coefficient 
drag 

= diameter 


Received June 23, 1958. 

+ This paper presents the results of one phase of research! 
carried out at the Jet Propulsion Laboratory, California Insti- 
tute of Technology, under joint sponsorship of the Department 
of the Army, Ordnance Corps (Contract No. DA-04-495-Ord 18), 
and the Department of the Air Force. 

** Chief, Aerodynamics and Propellants Department, Jet 


Propulsion Laboratory. 


Subscripts 


height of center of pitot tube above model surface 
root mean square average roughness height 
longitudinal length of cylindrical afterbody 
Mach Number 

nose 

number of sand grains per unit surface area 
Prandtl-Schlichting-Nikuradse data 

static pressure 

tunnel supply pressure 

total pressure measured by pitot tube 
dynamic pressure of fluid 

universal gas constant in equation of state 
Reynolds Number 

absolute temperature 

velocity in boundary layer in x direction 
velocity aty =k 

free-stream velocity 

‘friction”’ velocity 

coordinate along axis of model 

coordinate normal to axis of model 

ratio of specific heats (C,/C,) 
boundary-layer thickness 

boundary-layer momentum thickness 
viscosity 

kinematic viscosity 

fluid density 

shearing stress 


base 

critical 

fore 
incompressible flow 
skin friction 

total 

wall 

free stream 

local 


(1) INTRODUCTION 


HE EFFECT of surface roughness on the skin-friction 
drag of a surface is of great practical importance in 





bo 


the fields of aeronautical, marine, and hydraulic engi- 
neering. In the first two cases, it is related to the per- 
formance possibilities of aircraft and ships, and, in the 
third case, it is related to the efficiency of hydraulic 
machinery and installations. It is also of fundamental 
interest as a problem in fluid mechanics. 

Various aspects of the problem have been treated 
by many investigators in the low-speed incompressible 
case. The most systematic experimental investigation 
of both the skin-friction drag and the development of 
the boundary-layer structure on a rough surface was 
carried out by Nikuradse.” His tests were made with 
water flowing at low speeds through cylindrical pipes 
whose interior walls were artificially roughened by 
means of closely packed sand grains. More recently, 
experiments have been carried out at the Iowa Insti- 
tute of Hydraulic Research; this work has been re- 
ported in a very interesting paper by Hama.* The 
main results are concerned with the effect of rough- 
nesses of different kinds on the boundary-layer velocity 
profile. General skin-friction-drag formulas were de- 
duced from the velocity-profile laws, but no direct 
measurements of skin-friction drag were made. A 
summary of the present state of knowledge of the tur- 
bulent boundary layer for the low-speed case, including 
the effects of roughness, is to be found in the excellent 
paper by Clauser.* 

It is interesting to observe that, during the past 30 
years, research work on roughness effects has been 
piecemeal; the only really systematic investigation, that 
of Nikuradse,” was conducted with pipe flow. The 
resistance charts prepared by Prandt] and Schlichting? 
were developed from Nikuradse’s pipe data and involve 
the tacit assumption that the roughness effect would 
be universal and independent of outside-flow condi- 
tions. To this day, there do not exist (to the writer’s 
knowledge) any direct skin-friction-drag measurements 
for a low-speed boundary-layer flow on a rough sur- 
face. The writer has recently learned that a system- 
atic program of direct drag measurements on the 
flat wall of a low-speed wind tunnel with rough surfaces 
is now under way at the Aerodynamics Institute in 
Géttingen (1957). 

The one piece of work done at supersonic speeds on 
the drag of rough surfaces, to the writer’s knowledge, 
is that of Wade at the University of Toronto.’ Here, 
the critical roughness size which causes an increase in 
drag above that of a smooth surface was found for a 
screw-thread roughness at the single Mach Number 
of 2.48 and at one Reynolds Number. 

It was the purpose of the present investigation to 
determine experimentally the manner in which the 
turbulent skin-friction drag of rough surfaces varies 
with Mach Number and Reynolds Number in the 
supersonic-speed range and to compare these results 
with the available incompressible results. It was 
believed that such a comparison would show whether 
any new phenomena are involved in the compressible 
case. It was intended to find also the critical per- 
missible roughness of a surface up to which the skin- 
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friction drag is the same as that for a smooth surface. 
Finally, it was planned to make boundary-layer profile 
measurements in order to ascertain in what manner 
the boundary-layer structure is different from that in 
the incompressible case. 


(2) EXPERIMENTAL PROGRAM 


The experimental program consisted of two major 
parts. The first part was concerned with the direct 
measurement of skin-friction drag on rough surfaces. 
The second part was concerned with determining the 
effect of the surface roughness on the structure of the 


boundary layer flowing over the rough surface. 


(a) Skin-Friction-Drag Measurements 


(1) The JPL 18-In. by 20-In. Supersonic Wind Tun- 
nel—The skin-friction-drag tests were carried out in 
the 18-in. by 20-in. supersonic wind tunnel at the Jet 
Propulsion Laboratory. This tunnel is of the con- 
tinuous-flow type, has an 18-in. by 20-in. test section, 
and operates in the Mach Number range of about 1.3 
to 5.0. By varying the supply pressure, the Reynolds 
Number can be varied from about 120,000 to about 
400,000 per in. of characteristic length. The nozzle of 
the tunnel is formed by a pair of flexible steel plates, 
set to shape by 22 pairs of motor-operated jacks; thus, 
the tunnel may be run at any Mach Number in the 
range 1.3 to 5.0. The Mach Number is uniform in the 
test section to within +0.01, and the static pressure is 
uniform to within about +1 per cent. The flow is 
parallel throughout the test section to within about 
+(0.1°. At the present state of the art of wind- 
tunnel design, it can be said that this tunnel has a very 
high degree of flow uniformity. A photograph of the 
nozzle, looking upstream from the test section, 1s shown 
in Fig. 1. In this view, one of the side plates has been 
removed. 

(2) Measurement Technique—The 
for the drag measurements was selected because of its 
directness, simplicity, and adaptability to the wind 
tunnel and balance systems available. The models were 


technique used 


bodies of revolution, consisting of a smooth ogival nose 
piece followed by a cylindrical afterbody. The cylin- 
drical afterbody was carefully wrapped with a high- 
grade commercial sandpaper having sand grains of 
several different average sizes. A front-lighted schlie- 
ren photograph of a typical model, installed in the test 
section of the tunnel with the wind on, is shown in 
The drag-measuring technique was exactly 
The first step 


Fig. 2. 
that used by Chapman and Kester.’ 
was a test of the complete model, during which both the 
total drag D;, and the base pressure Dp, were measured 
(Fig. 3). From these two drag measurements, the 
fore drag of the complete model D;, was found: 


Dr, —_ Dz, == Dp, (1) 


Next, the smooth ogival nose was tested alone, the 
total drag and base pressure again being measured 
(Fig. 4). From these two drag measurements, the fore 
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The JPL 18-in. X 20-in. Supersonic Wind Tunnel. 





Fic. 2. Typical model installation. 





3. Schlieren photographs of complete model. 


(a) Left: 


M = 2.56. 


(b) Right: 


M = 0.70. 





Fic. 4. 
(a) Left: 


Schlieren photographs of ogival nose. 


M = 3.07. 


(b) Right: 


M = 0.70. 
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Fic. 5. Balance and suspension system of 18-in. X 20-in 
Supersonic Wind Tunnel (schematic). 





Fic. 6. Balance and suspension system of 18-in. X 20-in. 
Supersonic Wind Tunnel (photograph). 
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Fic. 7. NOL drag balance. 


drag of the ogival nose was then found: 
Dr, —_— Dz, = Dr, (2) 


Finally, by subtracting the fore drag of the ogival 
nose from the fore drag of the complete model, the re- 
sulting skin-friction drag of the cylindrical afterbody 
Dsr was found: 


Dr, — Dy, —_ Dsr (3) 


The average skin-friction-drag coefficient Cr of the 
cylindrical afterbody was determined, based on the 
wetted area of the cylindrical afterbody, and the corre- 
sponding Reynolds Number Ke was determined, based 
on the free-stream velocity V and the longitudinal 
length / of the cylindrical afterbody: 


Cr = (Dsr + AD)/[(p/2)V2ad (1 + Al)]) 
Re = V(l+ AD/» 


where p is the fluid density, d is the body diameter, and 
v represents kinematic viscosity. 

In all cases, the boundary layer of the ogival nose 
was made fully turbulent by means of a small wire 
ring placed at the nose of the ogive. In this fashion, 
it was made certain that the boundary layer over the 
rough cylinder was fully turbulent at all times. Both 
the trip wire near the nese and the resulting turbulent 
state of the boundary layer can be seen in Figs. 2, 3, 
and 4, 


> (4) 
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Because the turbulent boundary layer was thus of 
finite thickness at the front end of the rough-surfaced 
afterbody, it was necessary to find the effective length 
of the cylinder. This was accomplished by adding to 
the geometric length / a length A/ which was equal to 
that calculated fictitious ‘‘starting length’ of the cylin- 
der which would produce the actually observed bound- 
ary-layer thickness existing at the front end of the 
cylindrical afterbody. The drag corresponding to Al 
is represented by AD. In this way, both the average 
drag coefficient and the Reynolds Number were cor- 
rected to “‘effective’’ values, as shown in Eq. (4). 
The calculation of A/ and AD is given in Section (2.a.3). 

To form nine models, eight rough-surfaced cylindrical 
afterbodies and one smooth afterbody, together with a 
single smooth ogival nose piece, were used. For these 
nine models, the root mean square average roughness 
height k, in inches, was as follows: 


Q (smooth), 0.0008, 0.0027, 0.0039, 0.0066, 6.0095, 
0.0131, 0.0240, and 0.0380. 


Each cylindrical afterbody was undercut by an 
amount equal to the combined thicknesses of the sand, 
plus glue, plus paper (which combination makes up 
the sandpaper). Thus, the tops of the sand grains on 
the cylinder were a flush continuation of the smooth 
surface of the ogival nose at its base. 

Two balance systems were used during the investi- 
gation. The first of these was the standard external 
six-component mechanical-hydraulic balance with which 
This system, 


the wind tunnel is regularly equipped. 
which is built into the tunnel proper, is shown sche- 
A photograph of the balance out- 
The model is always 


matically in Fig. 5. 
side the tunnel is shown in Fig. 6. 





Fic. 8. 





Fic. 9. Sting-model connection and base-pressure annulus. 
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EFFECT OF UNIFORMLY 


at zero angle of attack and zero yaw; therefore, for the 
purposes of this investigation, only the drag component 
was used. The model drag pushes the entire balance 
and suspension system downstream, and this motion is 
resisted by the hydraulic load cell, marked D in Fig. 5. 
The model, sting, sector, sector bearings, moment 
table, pyramidal struts, force table, and drag link are 
all solidly connected and move downstream as a unit. 
The drag link connecting the force table to the drag cell 
D can be seen in Fig. 5. The performance of this bal- 
ance is very reliable and repeatable; its accuracy is 
discussed in Section (2.a.5). 

Since the skin-friction drag is found from a double 
subtraction of four measured quantities, and since there 
is a large change in dynamic pressure as the Mach 
Number is changed from 2 to 5 and the Reynolds 
Number is varied by a factor of 4, there was some con- 
cern at the start of the investigation about the accu- 
racy of a fixed-range balance, such as the JPL balance 
described above. Therefore, at the beginning of the 
program, an especially designed single-component vari- 
able-range spring balance* was used. Sketches of the 
Naval Ordnance Laboratory balance are shown in 
Fig. 7, and its operation is fully described by Kendall.* 
It will be seen that it is simply a direct spring balance. 

* This balance, which was designed and built by J. M. Kendall 
of the U.S. Naval Ordnance Laboratory (NOL) in White Oak, 
Md., was known to the writer because of its successful use at 
NOL. It was loaned to JPL for an extended period of time by 
Mr. Kendall and Dr. H. Kurzweg (of NOL). In addition, Mr. 
Kendall spent a week at JPL, during which time he instructed 


Laboratory personnel in the use of the balance. The author 
would like to express here his grateful appreciation to Dr. Kurz- 


weg and Mr. Kendall for their kindness and assistance. 
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Fic. 10. Schematic diagram of reader. 
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Fic. 11. 


Reader box and associated electrical equipment. 


The sliding friction is minimized by the use of an electric 
motor to rotate the bearings which support the model 
The balance load range is changed by inserting 
The axial de- 


sting. 
a spring of the appropriate stiffness. 
flection of the sting is measured by the stem attached 
to the sting frame. The stem is located inside the 
field of a Schaevitz coil, which electrically picks up its 
motion. A reader box outside the tunnel contains a 
second Schaevitz coil and a finely divided micrometer, 
which, through a lever and gears, moves a stem in the 
coil. With no drag applied, the two coils are elec- 
trically nulled. When the tunnel airflow is turned on 
and the drag is applied, the force extends the spring 
and upsets the balance between the two coils. The 
micrometer is then rotated manually until the outputs 
of the two coils are again nulled and the micrometet- 
drum rotation to balance is proportional to the drag. 
The entire system is mechanically and electrically 
calibrated. 

The model, with sting and balance box, is shown in 
Fig. 8. Fig. 9 shows a close-up view of the sting-model 
connection and the annular gap through which the base 
pressure is led inside the model; from here it proceeds 
through the sting to the balance box and to the balance- 
sector windshield, whence it goes outside the tunnel 
to a manometer for recording. 

Fig. 10 presents a schematic diagram of the reader. 
Fig. 11 shows the reader box, containing the micrometer 
and the second coil, and also shows the associated 
electrical equipment needed for power supply, and for 
voltage and voltage-balance indication between the two 
coils. 

(3) Data Reduction and Corrections 
steps in the data-reduction process have been indicated 
in Eqs. (1) to (4). For all models, / = 14 in. and 
d = 2in. The effective ‘‘starting-length’’ correction 
A/, due to the finite boundary-layer thickness at the 
front end of the cylindrical afterbody, and the skin- 
friction-drag correction AD, corresponding to Al, are 
found as follows. 

By equating the sum of the local skin-friction drag 
up to a station x on a surface to the loss of momentum 
in the boundary layer at x, the following expression is 
obtained : 


The essential 
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EFFECT OF UNIFORMLY 


f bz > 
D(x) = ad iy T(x) dx = mdp f, u(V — u) dy 


Here, x designates the coordinate along the axis of the 
model, y represents the coordinate normal to the axis 
of the model, 7 is the shearing stress, 6 is the boundary- 
layer thickness, and wu is the velocity in the boundary 
layer in the x direction. But the boundary-layer 
momentum thickness 6 is given by 


P= v) | u(V — u) dy 
( 


and, therefore, D(x) = mdpV0 


Now, Cr = D(x)/[(p/2) V2adx] = 

wd pV70/|(p/2) V*xdx|] = 20/x 
Therefore, Al = 20/Cr (5) 
and AD md pV"6 (G6) 


These are the corrections in Eq. (4). 
In the present investigation, the boundary-layer 


thickness § at the forward end of the cylindrical after- 
body was measured optically, and 6 was found from 


d= O measured (6 6) 


The denominator is the calculated ratio 6/6 for a 
1,/7-power-law turbulent boundary layer at supersonic 
speeds.” 

In Eq. (5), Cpr is the average turbulent skin-friction- 
drag coefficient for a rough surface at a Reynolds 
Number based on the length A/ and on the velocity at 
the edge of the boundary layer. 

For the models used, the length / was 14 in. The 
correction A/ in Eq. (4) varied from 3 to 9 in., depend- 
ing on the Mach Number, and AD was from 20 to 50 
per cent of Ds,y, also depending on the Mach Number. 

The present tests were made using ogive-cylinder 
bodies of revolution because of the intrinsic ease and 
reliability of testing with such models, rather than with 
a two-dimensional model such as a flat plate. For long 
cylinders, the boundary-layer thickness becomes com- 
parable to the cylinder radius; under these conditions, 
appreciable departures from flat-plate skin-friction 
values would be expected, according to the analyses 
of Jakob and Dow’ and Eckert.'' Chapman and 
Kester’ showed experimentally that the skin-friction- 
drag coefficient for a model with an //d value of 23 was 
a little less than 5 per cent higher than for a model 
with an //d value of 8. Furthermore, Chapman and 
Kester obtatmed the same skin-friction-drag results 
with bodies of revolution as did Coles,!* using a flat 
plate. (See Chapman and Kester,’ Fig. 11.) Finally, 
the //d value for the models of the present investigation 
was 11.5 in all cases; it is therefore concluded that the 
/‘d value for the present models is sufficiently small 
that the results are the same as those for a flat plate. 

The calculation of the roughness Reynolds Number, 
shown in Section (3) to be important, is as follows: 


Roughness Reynolds Number = kvx/y 


where vx = “‘friction’’ velocity = a/r p. 


DISTRIBUTED ROUGHNESS r 


In calculating kvy /v, the wall values of all factors are 
used, indicated by the subscript w; the subscript 0 
designates a free-stream quantity; and the local value 
of the skin-friction-drag coefficient is indicated by C,: 


/ 
Vx (Vv = kv Tw) Pw) Vu 
RV Cy(po/2) Vo?/ pw/¥m = RVo VC,s/2 V po/fw/ Ve 


Now, assuming that Prandtl Number = 1; p = pRT, 
where p static pressure, RK = universal gas constant 
in the equation of state, and 7 = absolute temperature; 
p is constant across the boundary layer; and u ~ TJ”, 


where uv = viscosity and w = 0.76; and setting 


1/}1 + [(y — 1)/2)M*?} = A 


where y = ratio of specific heats; then 


kuz /v = RVoVC,/2V1 


kve/v = RVoWVC,/2/(vf1 + [Cy — 1)/2]M?} 


Here, C, corresponds to a given k and M. 

(4) Skin-Friction-Drag Results—The results of the 
drag measurements described above are plotted in 
Figs. 12 to20. The data in Figs. 12 to 17 were obtained 
using the JPL balance, and the data in Figs. 18 to 20 
were obtained using the NOL balance. In each of these 
graphs, the average skin-friction-drag coefficient is 
plotted against the Reynolds Number, based on the 
effective length of the cylindrical afterbody, for several 
values of the roughness size at a given Mach Number. 
For purposes of comparison, each graph includes the 
skin-friction drag of an aerodynamically smooth surface 
for the same Mach Number, plotted according to the 
measurements made by Coles in the same wind tunnel 
in 1952.1 It is of immediate interest to compare the 
present measurements for a rough surface with the same 
measurements for a rough surface at / 0. It is 
possible to do this by using Nikuradse's data,? which 
are discussed in Section (1). These results were con- 
verted to the drag coefficient for a flat-plate boundary 
layer by Prandtl and Schlichting. They are shown 
in Figs. 12 to 17 and are identified by the symbol 
P-S-N (M = C). It can be seen immediately in Figs. 
12 to 20 that, at supersonic speeds, the drag coefficient 
decreases steadily as the Mach Number is increased. 
A comparison between the present supersonic results 
and Nikuradse’s results at 1/ = 0 is made at the same 
values of //kR. 

Fig. 21 shows the compressibility effect at supersonic 
speeds for rough surfaces. In this graph, the results 
shown in Figs. 12 to 17 have been crossplotted at a 
Reynolds Number of 8 X 10°. The ratio of the com- 
pressible skin-friction coefficient, as determined in the 
present investigation, to the incompressible skin- 
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friction coefficient, as determined by Nikuradse, is 
plotted against Mach Number; it is seen that the re- 
sulting curve is fair, and that the variation with Mach 
Number is similar to that determined experimentally 
by Coles!” for the case of smooth surfaces. It is also 
apparent that the compressibility effect for rough sur- 
faces is considerably greater than that for smooth sur- 
faces. Asa matter of interest, Fig. 21 also includes the 
estimate of the compressibility effect for smooth sur- 
faces made by von Karman in his 1935 Volta Congress 
paper.'* It is important to point out that the skin- 
friction-drag ratio at each of the six Mach Numbers 
of the present tests was independent of roughness size. 
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In Fig. 22, the compressibility effect is shown in a 
more directly physical way, which was suggested by 
Liepmann.'* If the roughness height is such that the 
skin-friction drag is indeed quadratic, then one may 
write 


Cr = Conk(1/2) pyu,?/[(1/2) pV?] 


where Cp is the drag coefficient of a single sand grain 
and is nearly independent of Reynolds Number and 
Mach Number, since the roughness element is a bluff 


body; m = number of sand grains per unit surface 
area; and u = velocity at y = k. 
Also, 


Cr, = Cpnk(1/2) pu,?/[(1/2)pV?| 


where Cp; is the average skin-friction-drag coefficient 
Then, to first order, 


2|M?| 


for incompressible flow. 


Cr Cr; = fel? = l {1 + 0.86[(7 = 


In this manner, the data of Fig. 21 have been plotted 
in Fig. 22 and compared with the straight-line relation- 
ship 


Cr Cr, = Pw’ P 


Indeed, the comparison is excellent. From this, one 
can see that the skin-friction drag of the rough surface 


is quadratic, and that the entire compressibility effect is 


a reduction of the density at the wall as the Mach 
Number increases. 
Thus, the density effect, which is largely offset by 


viscosity effects for laminar skin friction and partly 
offset by viscosity effects for turbulent skin friction on 
smooth surfaces, is completely dominant for rough sur- 
faces. Hence, the decrease friction with in- 
creasing Mach Number is most pronounced when the 
surface is fully rough." 

It would be interesting to see how the ratio of the 


in skin 


compressible skin-friction coefficient to the incompres- 
sible skin-friction coefficient varies as the roughness 
size is steadily decreased from the hydraulically smooth 
condition to the aerodynamically smooth condition. 
This is shown in Fig. 23, using the data of Fig. 24. 
The data in Figs. 23 and 24 are taken from measure- 
ments made by Wade’ at the single Mach Number of 
2.48 for very small roughnesses. It can be seen that 
the drag-coefficient ratio increases from the roughness 
size which is critical down to the aerodynamically 
smooth case, and that, when k = 0, the drag-coefficient 
ratio is the same as that determined by Coles for a 
smooth flat plate. 

The data of Figs. 23 and 24, 
lowing important observations: 

(a) Wade’s “‘roughness’”’ was a screw-thread rough- 
ness; that is, it was obtained by cutting a screw thread 
of peak-to-valley height k on the cylindrical portion 
of his model. Fig. 24 shows that the drag increase due 
to a sand-grain roughness (JPL data from Fig. 14) is 
about the same as that due to a transverse-screw- 
thread roughness of the same height k. The screw- 
thread roughness is two dimensional. 


in fact, lead to the fol- 
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(b) Wade’s model was a cone-cylinder model, and, 
for the smooth case (k = 0), he obtained the same com- 
pressible-to-incompressible skin-friction-drag ratio as 
did Coles for a flat plate (Fig. 23). Thus, an inde- 
pendent check is provided on the previous demonstra- 
tion that, for suitable length-diameter ratios of bodies 
of revolution, the skin-friction drag is the same as that 
on a flat plate. 

(c) Fig. 24 shows that, at WM =~ 2.5, the critical 
roughness at which the surface becomes hydraulically 
smooth is about k = 0.0008 in.; thus, the values of 
Cr/Cr, in Fig. 23 which are intermediate between 
k = 0 and the asymptotic value for large k correspond 
to hydraulically smooth surfaces in the compressible 
case. Finally, this means that the area between the 
“rough” and “‘smooth’’ compressibility-effect curves 
in Fig. 21 corresponds to surfaces which have become 
hydraulically smooth in the compressible case (Cp), 
but are still hydraulically rough in the incompressible 
case (C,r;) for the same value of //k. 

In Fig. 25, the skin-friction-drag coefficients from 
the present tests have been plotted against Mach 
Number at the Reynolds Number of 8 X 10° and are 
compared again with the values for a smooth surface 
as determined by Coles. It is clear that, as the Mach 
Number increases, a surface of a given roughness be- 
comes hydraulically smooth because of the growth of 
the boundary-layer thickness with increase in Mach 
Number. 

In Fig. 26, the skin-friction-drag coefficient is plotted 
against the roughness height at a Reynolds Number 
of S X 10° tor each of the six Mach Numbers of the 
present tests. 

Fig. 27 presents a correlation of the results of Fig. 
26; this is obtained by plotting the ratio of the skin- 
friction-drag coefficient for the rough surface to the 
coefficient for the smooth surface versus the logarithm 
of the roughness Reynolds Number kv / v, based on the 
roughness height k and the friction velocity v¥. The 
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obviously successful correlation shown in Fig. 27, 
Cr/Cr, = fllog (kv /v) |, is the same as the functional 
relationship shown to hold in the incompressible case 
by Nikuradse? for closely packed sand-grain roughness. 
(5) Accuracy of the Measurements—The accuracy 
problem, as well as the accuracy actually achieved in 
these present tests, is demonstrated by considering 
two extreme cases from the present investigation: 


Case (1) 


M = 1.98 
k = 0.038 in. 
Pr = 140 cm. Hg (supply pressure) 
Dy, = 0.061 X 140 = 8.53 lb. 
Dr, = 0.023 X* 140 = 3.22 lb. 
Dsr = 8.53 — 3.22 = 5.31 Ib. 
If p, = 45 cm. Hg (a minimum value at WM = 1.98), 
Case (2) 
M = 4,54 
k = 0.0095 in. 
pr = 300 cm. Hg 
Dr, = 0.0036 X 300 = 1.68 lb. 
Dr, = 0.0026 X 300 = 0.78 Ib. 
Dsr = 1.08 — 0.78 = 0.30 lb. 


If p, = 140 cm. Hg (a minimum value at M = 4.54), 
Dgpr = 0.16 lb. 


Now, the JPL external six-component balance has 
a low-scale drag range of +20 lb. and, at this range, 
has an accuracy of about +0.02 Ib. It is clear from 
the examples given above that, at low Mach Numbers 
and large k values, the accuracy is a fraction of 1 per 
cent, and that, with increasing Mach Number and de- 
creasing roughness, the accuracy becomes steadily 
worse until, in the worst case, it is about 10 per cent at 
M = 4.54andk—0O. For the majority of the data, it 
is estimated that the overall accuracy is about 2 to 5 
per cent. It is believed that this is adequate, both fot 
any practical engineering use of the results and for 
numerical proof of the physical-aerodynamical con- 
clusions which are drawn in Section (4). 

It is interesting to observe in Fig. 17 (JPL balance) 
and Fig. 20 (NOL balance) that the internal consistency 
of the data during a continuous wind-tunnel run is 
about an order of magnitude better than the accuracy 
from run to run; this would be expected. The k = 0 
values in Fig. 20 show the greatest scatter experienced 
from run to run under the most difficult conditions 
encountered—1.e., maximum M/ and minimum k. The 
mean curve drawn through the data for k = 0 is about 
7 per cent below the data of Coles, which are certainly 
the most accurate data for smooth surfaces available. 


(b) Boundary-Layer Measurements 


(1) The JPL 12-In. by 12-In. Supersonic Wind Tun- 
nel—Velocity-profile measurements were made in the 
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The JPL 12-in. X 12-in. Supersonic Wind Tunnel, 


Fic. 28. 


turbulent boundary layer of a cone-cylinder model, the 
cylindrical portion of which was wrapped with sand- 
paper having sand grains of several sizes. The tests 
were conducted in the 12-in. by 12-in. supersonic wind 
tunnel of the Jet Propulsion Laboratory. This wind 
tunnel is continuous, and has a Mach Number range of 
about 1.3 to 4.0; the Reynolds Number is also variable. 
The jack-operated flexible-plate nozzle and test sec- 
tion, with a model in place, is shown in Fig. 28. The 
flow in the test section is very uniform and is quanti- 
tatively the same as that given in Section (2.a) for 
the 18-in. by 20-in. tunnel. 

(2) Measurement Technique 
for the boundary-layer profile measurements is shown 
The cone-cylinder model 


The technique used 
in the photograph of Fig. 29. 
vas mounted rigidly in the center of the tunnel, parallel 
to the air stream. A small pitot tube was mounted on 
a vertical-traversing strut, which can be hand-cranked 
by means of a rack-and-pinion gear drive down through 
an opening in the ceiling of the test section. The in- 
ternal opening of the pitot tube was rectangular, with 
dimensions of 0.005 in. vertically and 0.012 in. later- 
ally. Outside dimensions were 0.008 in. and 0.018 in. 
The entrance of the pitot tube was fabricated from 
razor-blade pieces. 

A test run was made by cranking the pitot tube down 
in one direction (to remove the gear backlash) through 
the turbulent boundary layer until the tube touched 
the surface of the model. Contact with the model 
surface or with the top of the sand grains was observed 
optically through the tunnel window, using a Wild 
telescope level. This could be done with good accu- 
racy, since the flow, model, and probe were all solid and 
stationary. 

A typical schlieren photograph of the flow field and 
the boundary-layer development along the model is 
shown in Fig. 30. Here, again, a wire ring at the nose 
of the model produced a turbulent boundary layer 
along the cylindrical portion of the model. Weak Mach 
lines can be seen emanating from the distributed 
roughness starting just downstream of the shoulder. 

(3) Data Reduction—The measurements of u versus 
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the distance y from the body are put into a nondimen- 
sional relationship in order to compare them with the 
corresponding incompressible results. 

For flow near the wall, the relation is the Prandtl 
wall velocity law for incompressible flows, 


u/ve = f(yve/v, ¥/k), (y > 0) (8) 


Now, 
Vx /b = yV 1% Pu/Yo = 
yViV C,/2/(nfl + (Cy — 1)/2]M2} 2) 9) 
This equation was derived in Eq. (7). Also, 
u/V, = U V ty a, = (/ Vy) X 
C/V C,/2 {1 + [(y — 1)/2] M2} ) (10) 


Here, the subscript 1 indicates values at the outer edge 
of the boundary layer, where y = 6. 

The nond:mensional ratios in Prandtl’s relation above 
have thus been written in terms of the values of the 
variables at the wall for the compressible fluid. 

A generalization of the incompressible -mixing-length 
theory to compressible flow by van Driest'® shows that 
u/vx Should be replaced by the term 


sin~! m(u/1) 
MYV Tw/ Pw/ U1 


and that yu, /v should be replaced by 


/ 
WV Tw/ Pw/ Vy 
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Model and pitot-tube installation for boundary-layer 
measurements. 
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Schlieren photograph of cone-cylinder mcdel in 12-in. 
X 12-in. tunnel. 
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In the relation 


sin ~!m(u/u) y W tt Po Y 


“sy x 11) 
mV Pw/ Uy Vy k 
m*? = (y — 1)M,?/[2 + (v — 1)M?] 


Eqs. (9), (10), and (11) permit the complete data 
reduction of the results of the boundary-layer measure- 
ments discussed in following paragraphs. 

(4) Boundary-Layer Results—Boundary-layer meas- 
urements using the total-head tube technique described 
above were made on a smooth model (k = 0) and on 
models with sand-grain roughnesses of k = 0.0027, 
0.0052, and 0.0131 in. The measurements were made 
at a station 8 in. aft of the nose, as shown in Fig. 29, at 
a free-stream Mach Number .) = 2.60, with a corre- 
sponding local-surface Mach Number M, = 2.75. 

The raw data for the smooth model (k = 0) and for 
the three rough models are plotted in Figs. 31 to 34. 
The ratio of the total pressure measured by the pitot 
tube p,’ to the supply pressure p, is plotted against 
the height / of the center of the pitot tube above the 
surface of the model. 

The velocity profiles obtained assuming constant 
total energy and constant static pressure through the 
boundary layer are shown in Fig. 35 where the ratio 
of the velocity u at a point y in the boundary layer to 
the velocity u; at the edge of the layer (y = 5) is plotted 
against y/6. 

The nondimensional velocity profiles are plotted in 
Fig. 36 in the form 


sin—! m(u/u) VV tT. Pw 


‘ versus 
m V+. Pw/ M1 Vy 
and in Fig. 37 in the form 
Uu tf +. Py» Versus yV 1% Din! Bag (13) 


The smooth-model data fit the incompressible loga- 
rithmic law for smooth surfaces, u/vx = 5.5 + 5.75 log 
(yvx/v), better when put in the form of relation (12) 
(see Fig. 36) than when presented in the form of rela- 
tion (13) (see Fig. 37). 

It is seen that the effect of roughness is to displace 
the velocity-profile curve downward to a position essen- 
tially parallel to its original position by an amount 


(12) 


sin! m(u/1) 
m V Tw. Pw/ U1 
which is a function of the roughness Reynolds Number 
kvx Vv. 
When a surface is fully rough—that is, when the skin- 
friction drag is quadratic—then Prandtl’s wall law 


[Eq. (8) | simplifies to 


u/ve = fily/k) (yO) (14) 


In the incompressible case, Nikuradse? found this law 
to be 


u/vx = 5.75 log (y/k) + 8.5 (15) 
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The velocity profiles of the present investigation are 
plotted in Fig. 38 in the form 


sin—! m(u/1u) 


versus log y/k (16) 
m W t% Pw/ Uy 
and in Fig. 39 in the form 
u/V Tw/ pw versus log y/k (17) 


In both cases, Nikuradse’s relation of Eq. (15) is included 
for comparison. 

The downward displacement of the velocity profiles 
of Figs. 36 and 37 asa function of the roughness Reyn- 
olds Number kvx/v is plotted in Fig. 40. It is clear 
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that the quantitative shift in the velocity profile, 
whether defined as in Fig. 36 or as in Fig. 37, is re- 
markably close to the value obtained by Nikuradse in 
the incompressible case.” 


(3) DISCUSSION AND INTERPRETATION OF RESULTS 


(a) Skin-Friction-Drag Results 


The reliability of the experimental technique used 
in this investigation is demonstrated in Figs. 18 to 20 for 
supersonic speeds. It is shown in these graphs that 
the measurements made with a smooth model are 
within 10 per cent of those obtained by Coles!* by 
integrating carefully made local skin-friction measure- 
ments to obtain the average skin-friction coefficients. 
In Section (2.b.5), it is shown that the combination of a 
smooth model and a high Mach Number produces the 
smallest drag forces and the largest errors. Thus, it 
may be concluded that the measurements made with 
rough models at k > 0.008 in. are accurate to within 
about 5 per cent, at least. 

At the single subsonic Mach Number (7 = 0.70) at 
which tests were made, the smooth-model data coin- 
cide with the skin-friction laws of von Karman" and 
Coles” for M = 0 within the experimental error of the 
present test data (Fig. 12). Now, the measurements 
of Dhawan! show that, at 1/7 = 0.70, the skin-friction- 
drag coefficient is about 3 per cent smaller than that 
of von Karman for M@ = 0. Thus, the present meas- 
urements at 7 = 0.70 are consistent with those of 
Dhawan at the same Mach Number. 

A comparison of the present data at = 0.70 with 
those of Nikuradse at WM = O (Fig. 12) shows that, 
for k = 0.024 in. and k = 0.038 in., the skin-friction- 
drag coefficient is, in each case, equal to about 89 per 
cent of Nikuradse’s values at 1J = (0. This com- 
pressibility effect is consistent with the supersonic 
compressibility effect for rough and smooth surfaces 
shown in Fig. 21. This general consistency between 
the present tests at 1J = 0.70, those of Dhawan at 
M = 0.70 for smooth surfaces, and those of Nikuradse 
at M = 0 for rough surfaces tends, in a general way, to 
confirm the correctness of the Prandtl-Schlichting- 
Nikuradse charts.’ The present data at WM = 0.70 for 


k = 0.0095 in. coincide exactly with the Prandtl- 
Schlichting-Nikuradse data at M = 0, as shown in 
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Fig. 12 for the same value of //k. Because all the 
JPL data correlate successfully at all values of k (Fig. 
27), it is concluded that the Prandtl-Schlichting- 
Nikuradse data for //k = 1,680 are about 10 per cent 
too low. 

The compressibility effect for rough surfaces, shown 
in Fig. 21, is seen to be considerably greater than for 
smooth surfaces and to be essentially independent of 
the roughness size k. 

The results of the present investigation of the aver- 
age skin-friction-drag coefficient are summarized in 
Fig. 26 at a Reynolds Number of 8 X 10°, which is 
the maximum value attained. 

In Fig. 27, the data of Fig. 26 are plotted as the 
ratio of the skin-friction-drag coefficient for the rough 
surface Cyr to the value for the smooth surface Cr, 
versus the logarithm of the roughness Reynolds Num- 
ber kvy/v. The correlation is seen to hold for al] 
Mach Numbers up to about 5 and is the same func- 
tional relationship as that shown to hold in the incom- 
pressible case by Nikuradse.” The important result 
follows that the Mach Number is eliminated as a vari- 
able per se; hence, the effect of surface roughness on 
skin-friction drag is localized deep within the boundary 
layer at the surface itself and is independent of the 
external flow. Now, this result (the roughness effect is 
universal and independent of the outside flow field) is 
precisely the assumption made by Prandtl and Schlicht- 
ing’ in preparing their resistance charts for flat-plate 
boundary-layer flow from Nikuradse’s data obtained 
from the flow of water inside circular pipes.” The 
present results obviously strengthen this assumption, 
as does the consistency of the subsonic data of the 
present investigation with that of Prandtl-Schlichting 
in Figs. 12 and 21, and also the consistency of the sub- 
sonic and supersonic data of Fig. 27. 


(b) Critical Roughness Size 


A further important result is obtained from the cor- 
relation of Fig. 27: the critical roughness Reynolds 
Number, below which the skin-friction drag of a rough 
surface is the same as that of a smooth surface, is k,vx /v 
~ 10 for all values of J up to about 5. Here, k, is the 
Hence, if k < 10y/vx, then the 
The practical im- 


critical roughness size. 
surface is hydraulically smooth. 
portance of this result is obvious. 

It is interesting to note that the value given above 
for the critical roughness size (k, ~ 10v/vx) is numeri- 
cally the same as the thickness of the laminar sublayer 
of the turbulent boundary layer on a smooth surface 
at both subsonic and supersonic speeds. 

Now, since the surface is hydraulically smooth if 
k is less than the laminar-sublayer height of a physically 
smooth surface, one is led back to the simple idea, 
suggested many years ago, that the mechanism of the 
drag increase due to roughness of the sand-grain type 
below the critical-roughness size the flow 
As the rough- 


is as follows: 
about the roughness element is laminar. 
ness size increases, a point is reached at which vortices 
are formed behind the roughness element; turbulent 
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EFFECT OF UNIFORMLY 
flow about the roughness element ensues, and the drag 
rises steadily acccrding to the quadratic law as the 
roughness size is increased further. Clearly, then, one 
is interested in the critical Reynolds Number, based on 
the roughness size and on the local velocity at the top 
of the roughness element. 
It is easy to show that, if k, ~ 10v/v%, then 


ku. /v = (10/VC,/2) (u/V) 


Since k, is about equal to the laminar-sublayer height 
for a smooth surface, u; is the velocity at the edge of the 
laminar sublayer. 

In order to make an estimate of k,u,/v, the quantity 
wu, V for a turbulent boundary layer at an average 
Mach Number of 3 is taken as about 0.65, and C; is 
taken as about 0.0015 at Re = 8 X 10%. Therefore, 


Ruy vy ~~ 250 


In a very interesting investigation of the effect of 
sandpaper roughness on boundary-layer transition at 
low speeds, Smith and Clutter'’ found that turbulent 
spots occurred in the laminar flow at the roughness 
when ku;,/v = 275; furthermore, this value was found 
to be critical. That is, if ku,/» was slightly less than 
275, the flow in the boundary layer was completely 
laminar. The extent of the area covered by sand 
grains in the stream direction did not have any im- 
portant effect on the critical value of the roughness 
Reynolds Number. 

From the discussion presented above, it appears 
reasonable to conclude that the mechanism of the drag 
rise due to sand-grain roughness in turbulent flow at 
the critical roughness Reynolds Number and _ the 
mechanism of the critical occurrence of transition in the 
laminar boundary layer on a_ sand-grain-roughened 
surface are probably the same, and also that the 
mechanism is simply the breakdown of local laminar flow 
about the sand-grain-roughness element to form turbt- 
lent flow at a critical roughness Reynolds Number. 


(c) Velocity Profiles 

In regard to the velocity-profile results, it remains to 
observe that the turbulent-velocity profile for the 
smooth surface at supersonic speeds follows the Prandtl 


wall law 
uve = fllog(yx /v) | 


in its logarithmic portion (Fig. 36). 

When the surface is fully rough (that is, when the 
skin-friction drag is quadratic), the velocity profiles 
follow the Prandtl wall law 


u/vx = fllog(y/k) | 


in their logarithmic portions (Fig. 38). 

The displacement A(u/v*) of the velocity profiles 
of Fig. 36 with increasing roughness size k is shown in 
Fig. 40 to be a function of only the roughness Reynolds 
Number kvx/v. It is clear, furthermore, that the 
present data at \/, = 2.75 check exactly with the low- 


speed data of Nikuradse.* Since the displacement 
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A(u/v) of the velocity profile is proportional to the 
drag increase due to roughness, and since the drag 
increase due to roughness is a function of only the rough- 
ness Reynolds Number (Fig. 27), it is clear that the 
exact check between the compressible and incompres- 
sible velocity-profile shifts shown in Fig. 40 must 
follow. Finally, then, the drag correlation of Fig. 27 
and the velocity-profile correlation of Fig. 40, taken 
together, completely unify the subsonic and super- 
sonic skin-friction-drag measurements of the present 
investigation and the incompressible skin-friction-drag 
measurements of Nikuradse.* ® 


(4) CONCLUSIONS 


The following conclusions can be drawn from the pres- 
ent investigation : 

(1) At subsonic and supersonic speeds up to a Mach 
Number of 5, and for roughness sizes such that the 
quadratic resistance law holds, the compressibility 
effect is indirect, and the skin-friction drag is a function 
of only the roughness Reynolds Number kv /v, exactly 
as in the incompressible case. 

(2) The entire compressibility effect is a reduction of 
the fluid density at the surface as the Mach Number 
increases. Thus, the density effect, which is largely 
offset by viscosity effects for laminar skin friction and 
partly offset by viscosity effects for turbulent skir fric- 
tion on smooth surfaces, is completely dominant for 
rough surfaces. Hence, the decrease in skin friction 
with increasing Mach Number is most pronounced 
when the surface is fully rough. 

(3) The critical roughness, below which the surface 
is hydraulically smooth, is k, ~ 10v/vs%, independent 
of Mach Number, and this is equal to the thickness 
of the laminar sublayer for a smooth surface for both 
compressible and incompressible flows. 

(4) Over the range of roughness sizes considered 
here, there appears to be no wave drag associated with 
the drag due to roughness. 

(5) The shift in the turbulent velocity profile 
A(u/vx) for a rough surface at supersonic speeds is a 
function of only the roughness Reynolds Number 
kvx/v and quantitatively follows identically the same 
law as that for the incompressible case. 
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On the Structure of Jets From Highly 
Underexpanded Nozzles Into Still Air 


T. C. ADAMSON, JR.,* anp J. A. NICHOLLS** 
The Unwersity of Michigan 


SUMMARY 


A method for calculating the position of the first normal shock, 
or Mach disc, in the jet behind a highly underexpanded nozzle is 
presented. In the calculation for a sonic orifice, the axial pres- 
sure distribution on the centerline of the flow behind the orifice, 
calculated by the method of characteristics, is used to define a 
fictitious nozzle extension, and the shock is then assumed to exist 
at that point where atmospheric pressure would be attained 
behind the shock—i.e., the shock is assumed to exist at the end 
of the fictitious nozzle extension. Physical arguments are then 
employed to extend this calculation to nozzles with supersonic 
exit Mach Numbers. The results compare favorably with experi- 
mental data. 

An approximate method for computing the jet boundary up 
to the point of maximum jet area is given, and the results are 
compared both with photographs of actual jets and with jet 
boundaries calculated by the method of characteristics. Favor- 
able agreement exists at relatively low nozzle pressure ratios. 


SYMBOLS 


A = area 

d = diameter of jet 

M = Mach Number 

P = pressure 

p = dimensionless static-to-total-pressure ratio, P/P 

r = radius of jet 

x = axial coordinate 

a = total divergence angle of jet boundary at nozzle exit 

6 = divergence half angle of nozzle 

& = dimensionless axial coordinate, x/dy 

y = angle made by jet boundary with axis of flow 

@ = total angle between the tangent to the jet boundary at 
any point and the original jet-flow direction, v — ve 

uw = Machangle, sin~! (1/1) 

vy = Prandtl-Meyer turning angle 

y = ratio of specific heats 

Subscripts 

t = total or stagnation conditions 

© = conditions in atmosphere into which nozzle is exhausting 

N = nozzle exit conditions 

e = conditions at jet boundary immediately after expansion 
to atmospheric pressure 

m = conditions at Mach disc location 


INTRODUCTION 


p | ‘HE INCREASED USE of jet and rocket engines and 
recoilless rifles has caused a great deal of attention 
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to be focused on the jet flow behind sonic and super- 
In particular, the flow behind highly 
A recent 


sonic nozzles. 
underexpanded nozzles is of great interest. 
research program conducted by Love and Grigsby at 
the NACA Langley Aeronautical Laboratory,' for ex- 
ample, includes comprehensive experimental and theo- 
retical studies of the structure of the jet flow behind 
nozzles of various configurations. 

Fig. 1 is a sketch of the repeating flow picture gener- 
ally observed behind an underexpanded axisymmetric 
nozzle. As the gas emerges from the exit, it goes 
through an expansion fan, expanding to the ambient 
pressure at the jet boundary, and then the condition 
of constant pressure along the jet boundary causes this 
boundary to be bent back toward the axis of flow. As 
the flow changes direction along this boundary, many 
compression waves, formed at the intersection of ex- 
pansion waves with the jet boundary, are sent back into 
the flow; these waves coalesce to form the intercepting 
shock indicated in the sketch. For slightly under- 
expanded nozzles, these intercepting shocks meet at the 
axis forming the familiar diamond configuration. As 
the pressure ratio across the nozzle is increased, how- 
ever, these intercepting shocks no longer meet at the axis 
but are connected with a normal shock, or Mach disc, as 
pictured in Fig. 1. In both cases, reflected shocks are 
formed which intersect with the jet boundary reflecting 
as expansion waves, and the whole process is repeated. 
The repetition is continued until viscous effects become 
predominant and this structure is no longer observed. 

The study of the underexpanded jet, then, involves 
calculations of the jet boundary, the intercepting shock, 
the position of the first normal shock, and the wave- 
length of the repeating structure. Experiment has 
indicated that for a jet exhausting into a static amos. 
phere the important parameters are the static pressure 
ratio at the nozzle exit (Py/P..), the Mach Number at 
the nozzle exit, and to a very slight degree, the diver- 
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gence angle of the nozzle. The pressure ratio, in par- 
ticular, has been used in several empirical relations 
for the wavelength and the distance between the nozzle 
exit and the Mach disc (see reference 1, p. 10, and 
reference 2). In this paper, only the downstream dis- 
tances to the normal shock and the jet boundary are 
discussed. Furthermore, simple approximate methods 
of analysis are presented which allow these calculations 
to be made quite rapidly. 


LOCATION OF THE First Macu Disc 


The downstream location of the first normal shock, 
or Mach disc, depends essentially on the static pressure 
ratio and the exit Mach Number, changing only slightly 
with the divergence angle of the nozzle for total diver- 
gence angles as high as 40°.'* This indicates that 
the viscous effects must be rather insignificant for this 
calculation; otherwise, the angle at which the flow enters 
the atmosphere would be an important parameter since 
the overall mixing area up to any downstream location 
increases with the initial boundary angle. Assuming 
viscous effects to be negligible, then, and noting that 
the main effects are due to the dynamic parameters 
with no geometric complications, one is led to the con- 
clusion that a simple quasi-one-dimensional calcula- 
tion should give accurate results. In this event, the 
problem is similar to the calculation of the shock posi- 
tion in a given nozzle, where isentropic quasi-one- 
dimensional flow is assumed up to and beyond the shock. 
This idea is also strengthened by a consideration of the 
picture usually observed when a shock is located 
within a nozzle. The normal shock does not extend 
from surface to surface, but extends only to an inter- 
section where a lambda configuration results; also, the 
flow separates from the wall slightly. However, in spite 
of this fact, the location of the shock is fairly well pre- 
dicted by disregarding the lambda configuration and as- 
suming the shock to be simply normal, extending across 
the whole fluid. Since the shock configuration occur- 
ring in the underexpanded jet is similar, although 
stretched out, it seems logical that there is a possibility 
of predicting the location of the normal shock without 
the necessity of calculating the whole shock structure. 
For the above reasons, then, the free jet was assumed 
to be simply an extension of the actual nozzle, bounded 
by a fictitious surface with an area distribution corre- 
sponding to the actual pressure distribution of the cen- 
ter, or core flow, up to the normal shock. The fictitious 
boundary extends only to the axial position of the 
shock wave, so as the total-to-atmospheric-pressure 
ratio changes, the length of the extension changes. 
Finally, since the pressure to which this fictitious nozzle 
is exhausting is atmospheric, the axial position of the 
shock is that point at which a shock has sufficient 
strength to bring the static pressure to atmospheric. 
In summary, the flow is assumed to be that through a 
It is assumed to ex- 
The nozzle 


nozzle with a shock at the exit. 
haust to atmospheric pressure at the exit. 
configuration is that which corresponds to a given axial 
pressure distribution in the center or core of the actual 


flow. Of course, the above description is simply a 
way of picturing the flow; it is not necessary to calculate 
the fictitious boundary to locate the shock since the 
axial pressure distribution along the centerline can be 
computed. The method of calculating the shock posi- 


tion is outlined in the following sections. 


Sonic Orifice 


Consider an underexpanded flow issuing from a con- 
vergent nozzle—i.e., Wy = 1. The problem of exit 
Mach Numbers other than one can be solved by a 
generalization of this case. The pressure variation 
along the centerline of such a flow has been calculated 
by Owen and Thornhill‘ for a total-pressure-to-atmos- 
pheric-pressure ratio, P,/P.., of infinity. However, as 
pointed out by Owen and Thornhill, the pressure vari- 
ation thus calculated should hold for all P,/P. ratios 
“in that region bounded by the orifice and the first 
wave front which registers the existence of an external 
pressure outside the jet.’’ Thus, for a considerable 
distance downstream, the pressure variation P/Py, 
where P is the static pressure at a given point along 
the centerline and Py is the static pressure at the 
nozzle exit, should be the same for all nozzle pressure 
ratios, P,/P 
region of similarity increases in size. 
hill plotted P Py vs. x/dy, where x is the distance 
downstream from the nozzle exit, and dy is the exit 
diameter of the sonic nozzle (Fig. 2). They also plotted 
the Mach Number at the centerline vs. x/dy (Fig. 3). 
Using these graphs, a curve of P2/Py vs. x/dy was 
plotted where P2 is the pressure behind a shock, occur- 
ring at a Mach Number corresponding to the given x /d 
(Fig. 2). It should be noted that the P/Py curve 
calculated by Owen and Thornhill was continued past 
their maximum x/dy of 3.85 so that Ps/Py could be 
calculated. Since the curve of M vs. x/dy (Fig. 3) 
was given up to x/dy = 10, this could be done simply 
by finding the P ‘Py corresponding to the Mach Num- 
However, since the P/ Py curve 


, and, as this pressure ratio increases, the 
Owen and Thorn- 


ber at a given x/dy. 
is so close to a straight line for x/dy > 1.2, and in view 
of the inaccuracies inherent in reading the original 
curves, a straight line extension was used, checking very 
well with calculations made from the Mach Number 
curve. 

The pressure behind the shock, P2, was taken to be 
Fig. 4 shows the results of such 
This cor- 


atmospheric (P..). 
calculations for a range of 2 < Py/P.. < 70. 
responds to a range of 3.8 < P,/P. < 133 and a max 
imum total pressure of approximately 2,000 psi if P. 
is standard sea-level pressure. On the same graph, the 
experimental results obtained both at the NACA 
Langley Field Laboratories and at the Aircraft Pro- 
pulsion Laboratory of The University of Michigan are 
plotted.* The error in x/dy between the theoretical 
and the NACA experimental curves is of the order of 
12 per cent to 14 per cent, while the error between the 

* The circles on the A.P.L. curve are actual experimental 
points; the NACA results are taken from curves drawn through 
their experimental points. 
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Fic. 5. Sketch of expansion waves emanating from the corner of 





an underexpanded nozzle with Wy > 1. 


NOZZLE “ 


Fic. 6. Sketch of expansion waves emanating from the corner of | 
an underexpanded sonic orifice 
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ON THE 
A.P.L. experimental and the theoretical curve is less. 
[he experimental setup for both laboratories is de- 
scribed in other reports,» * so they will not be repeated 
in this paper, except to note that, at the NACA, Love 
and Grigsby took schlieren photographs while at the 
A.P.L., shadowgraphs were taken. This could cer- 
tainly account for the difference between the two ex- 
perimental curves. Furthermore, the absolute error 
between experimental and theoretical values of x, the 
distance from the nozzle to the Mach disc, is of the 
order 0.07 in., which is certainly of the order of the ex- 
perimental error, the distances being measured on 
photographs. Therefore, the agreement between the 
theoretical and experimental curves appears to be most 
satisfactory, for the Mach one nozzle. 


Supersonic Nozzles 

For convergent-divergent nozzles with exit Mach 
Numbers different from one, it is necessary to find the 
pressure distribution, if the same method is to be used. 
To visualize better the following arguments concerning 
the pressure distributions, a sketch of the flow out of an 
axisymmetric convergent-divergent nozzle is given in 
Fig. 5. For simplicity in this figure and succeeding 
figures, only half the flow field is shown, and none of the 
intersection waves is drawn, although their effect is pic- 
tured through the curvature of the waves shown. The 
lines a—b, a—c, etc., are the wave fronts emanating from 


a~—b, in particular, 
1 


the expansion around the nozzle lip. 
is the Mach line defined by the wave angle, u = sin 
1 vy, i.e., it is the first wave in the expansion fan. It 
is apparent, from this picture, that, along the center- 
line, the fluid can be outside the nozzle, and still be at 
|My, the exit Mach Number of the nozzle. That is, 
until the centerline fluid reaches 5, no expansion signals 
However, after 6, the centerline receives 
which correspond to various Mach 
Again, the 


reach it. 
signals from a, 
Numbers at a, as the flow turns the corner. 
structure along the centerline should be the same for all 
total-pressure-to-atmospheric-pressure ratios until the 
wave front signaling the exterior pressure reaches the 
centerline. That is, the same argument holds for the 
centerline flow downstream from point b for the nozzle 
where Jy > 1 as holds for the centerline flow down- 
stream from the nozzle edge for the nozzle where My = 
1. However, as .\J/y changes, the distance from the 
nozzle to b changes and must be accounted for. Thus 
the problem is essentially that of finding the pressure 
distribution downstream from point 5. However, it 
is hypothesized that the desired pressure distribution 
can be found from the pressure distribution behind an 
My = The arguments sup- 
porting this hypothesis are as follows. 
Consider again the flow behind an J/y = 
As the expansion waves from the corner at 


1 nozzle (sonic orifice). 


1 nozzle 
(Fig. 6). 
a, intersect the centerline, the flow is expanded so that 
a variation of Mach Number exists. In fact, at a 
point b,, say, the Mach Number is the same as the Mach 
Number in the flow behind an J/y > 1 nozzle, at point 
b. Furthermore, at point a), as the flow turns through 
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an angle greater than that corresponding to the expan- 
sion wave d;—);, it expands through the same Mach 
Numbers, and thus passes through the same expansion 
waves as does the flow in Fig. 5, expanding beyond the 
line a—b. These waves propagate to the centerline in 
both cases, causing a variation in pressure and Mach 
Number, etc. This means that, as shown in Figs. 5 
and 6, there is a region on the centerline of both flows 
(b;-d, and b-d, say) where the same expansion takes 
place. The lengths },-d, and b-d are very nearly the 
same, unless the Mach Number of the nozzle, Wy, 
in the case where /y is different from one, is very high, 
when the initial angles of inclination of the lines a—)b and 
a,-b;, w and w, respectively, are very small. Then a 
small difference in starting angle, at a, and a, for any 
later expansion wave, means a large difference in dis- 
tance traversed on the centerline. The starting angles 
at a and a, for any of the expansion waves are differ- 
ent, of course, because in the /y > 1 nozzle the flow 
corresponding to the line a—d is parallel to the axis, while 
in the M/y = 1 case, the flow corresponding to line 
a,—b; is at an angle to the axis, depending on the value 
of the Mach Number associated with this line. 

The result of the above arguments is, then, that it is 
assumed that, for \/y > 1, the pressure distribution for 
x > x, is the same as the pressure distribution for the 
My = 1 nozzle flow for x > x,,. For a flow behind a 
nozzle with exit Mach Number of 1/y = 2, for example, 
b, is located at that point where the centerline Mach 
Number is 2. Furthermore, since the assumption is 
again made that the shock occurs so that the static 
pressure behind it is atmospheric, the distance from }, 
to the Mach disc in /y = 1 flow is the same as the dis- 
tance from } to the Mach disc in Jy = 2 flow, for the 
same total-to-atmospheric-pressure ratio, P,P... 

To describe analytically the above arguments, the 
following dimensionless notation is used: 


£) = (%n/dn) mya, &y, (Xp, dy) my =1 } 
Ag, = [(%m — Xs) dy lary ms £= dy)m \ (1) 
£, = (x,/d vV) ayy Aé = [ [ie = Ma dy | My 


where x,, is the distance from the nozzle to the Mach 
disc, and x, and x», are the distances between the nozzle 
and the points > and ,, respectively. dy is the exit 
diameter of either the sonic or supersonic nozzle de- 
pending on whether it is used in & or &. Point d is the 
point in the supersonic nozzle flow where the first wave 
of the expansion fan intersects the centerline, and point 
b, is the point on the sonic nozzle flow centerline where 
the Mach Number reaches the value of the exit Mach 
Number of the supersonic nozzle in question. For any 
given P,/P., & is a known quantity, found by the use 
of Fig. 2. Furthermore, for a given supersonic nozzle, 
i.e., for a given My, &, can be found from Fig. 3. Fi 
nally, knowing & and &, one can find Aé), since 


Afi = & — &, (2 


Next, & depends on yw, the Mach angle associated 


with J/y. 
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Fic. 8. Total angle turned through by the flow at the lip of an 
underexpanded nozzle, in terms of nozzle half angle and Prandtl- 
Mever turning angles. 


f= 12 tan ei (3) 


The value of x,, x, is the same as xy X»,, but to 
compute Aé,, Ag,, must be corrected to the exit diameter 


equivalent to the J/y of the supersonic nozzle. Thus, 


Ag, = Ab: (dn) ay =1/ (dn) my | ( 


Ag, = Ak, [(A*/Ae) my]? f (4) 


Then, one can write the equation for & as 
1/2 ‘ : 
& = Ag, [((A*/Ae)uy] ’~ + (1/2 tan p) (5) 


where Aé,, and & are calculated at the same P, P 

Eq. (5) was used to calculate the curves shown in 
Figs. 7(a), (b), (c), and (d), and the experimental values 
were also plotted for comparison. As in the case of 
My = 1, the comparison with experiment is quite good. 
Of course, the comparison between theory and experi- 
ment is made better by the fact that at low exit Mach 
Numbers the theoretical curve predicts values which 
are larger than experimental, while as the Mach Number 
increases, the difference becomes less at first and then 
increases so that the theoretical curve predicts values 
less than experimental curves. If the experimental 
values at M/y = 1 were used for &, the error at higher 
Mach Numbers would be considerably greater. 

It should be mentioned that several tests were car- 
ried out with various nozzle divergence half angles 
(6y), but the effect of nozzle divergence on the position 
of the first Mach dise was negligible, except perhaps at 
the higher nozzle Mach Numbers where the effects were 
still relatively quite small. 

It should be pointed out that if the dimensionless 
primary wavelength (the distance from the nozzle to 
the first constriction in the jet boundary, divided by the 
nozzle diameter—see Fig. 7, reference 1) is plotted ver- 
sus the static pressure ratio, Py/P.., on Figs. (7) of this 
paper, the slope of the curve is essentially the slope of 
the theoretical x,,/dy curve for each Mach Number. 
This indicates that the primary wavelength is a power 
function of the static pressure ratio, with the same 
power as the Mach disc position function, at each nozzle 
exit Mach Number. That is, the distance to the normal 
shock is proportional to the primary wavelength. How- 
ever, it appears that no further inference can be made, 
since the present equivalent nozzle hypothesis does not 
cover the flow past the normal shock. Furthermore, 
the primary wavelength seems to be much more af- 
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fected by the nozzle divergence angle than does the 
position of the Mach disc. 


CALCULATION OF THE JET BOUNDARY 


The ideal fluid jet boundary can be calculated using 
characteristic theory." °® This gives a good approxima- 
tion to the real jet boundary, although there is actually 
a viscous mixing region along the boundary. However, 
the work involved for each nozzle and pressure ratio is 
considerable and is still only approximate for the above 
reason and because there is some question about con- 
tinuation of the characteristic net beyond the inter- 
cepting shock. Therefore, the following method is 
presented as a relatively quickly calculated approxi- 
mation to the initial jet boundary. It does not agree 
with experimental photographs as well as the solutions 
obtained by the method of characteristics. However, 
the difference between the two becomes significant only 
at the higher jet expansion angles. At relatively low 
expansion angles, the agreement is quite good. 

Consider the flow at the lip of a nozzle (Fig. 8) with 
half angle dy. At the nozzle lip, before expansion, the 
Mach Number is \/y and the corresponding Prandtl- 
Meyer angle is vy. After expanding to atmospheric 
pressure, the Mach Number is \/, with a corresponding 
Prandtl-Meyer angle of v,. Thus, the flow at the nozzle 
lip turns through an angle of v, — vy relative to the 
nozzle wall, and the overall flow expansion angle with 
respect to the flow axis is a, where 


a= ve — vy + bn (6) 


The flow is turned from this initial expansion angle by 
the intersection of expansion waves with the boundary, 
reflecting as compression waves, so that the condition 
of constant pressure on the boundary is satisfied at 
every point. However, the present approximation 
depends on replacing the effects of expansion waves 
intersecting the flow near the boundary by the effects 
That is, 
just as in one-dimensional nozzle flow, it is assumed 
that one can find the average Mach Number and pres- 
sure at any axial position in an expanding flow from the 
area ratio at the given point, instead of going through a 
characteristic calculation. 

Immediately after leaving the nozzle exit, then, the 
flow at the lip has turned through the total angle a. 
However, as it follows this new direction, it is contin- 
ually expanding. Hence, in an incremental distance 
downstream, dx, there is an increase in area, dA, with 
a corresponding decrease in pressure, OP/O0A-dA. 
OP /OA can be calculated from quasi-one-dimensional 
relations if one considers a stream tube along the 
boundary of the jet to be the channel in question. 
Since a change in pressure along the jet boundary can- 
not be tolerated, this decrease in pressure must be bal- 
anced by an equivalent increase in pressure which can 
only be gained by turning the flow through an incre- 
mental angle, d@, thus forming a weak compression 
wave. This process of expansion and compression 
holds at any point on the boundary where at x, with a 


of a quasi-one-dimensional area increase. 
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Fic.9. Incremental angle and radius changes in the jet boundary 
due to an incremental increase in x, the axial coordinate. 
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Fic. 10(a). Comparison of jet boundaries calculated by ap- 
proximate method with those calculated by the method of char- 
acteristics for sonic and supersonic nozzles with various pressure 





ratios. (Characteristic solution after Love and Grigsby.') 
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Fic. 10(b). Comparison of jet boundaries calculated by ap- 
proximate method with those calculated by the method of char- 
acteristics for sonic and supersonic nozzles with various pressure 





ratios. (Characteristic solution after Love and Grigsby.') 
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Fic. 10(c). Comparison of jet boundaries calculated by ap- 
proximate method with those calculated by the method of char- 
acteristics for sonic and supersonic nozzles with various pressure 
ratios (Characteristic solution after Love and Grigsby.') 
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given 7, and a given direction of flow, @ (Fig. 9), the gas 
in expanding in this given direction to r + dr at x + dx 
would undergo a decrease in pressure so that a further 
change in direction, d0, results, keeping the pressure at 
a constant level at r + dr. 

The equation for the pressure along the jet boundary 
may then be written as follows: 

dP = 0 = (OP/0A)dA + (OP/06)d0 (7) 

Since a change in area, dA, is equivalent to a change 
in Mach Number, dV, Eq. (7) can be written in terms 
of /andé@. Thus, 

(OP/OM)dM + (OP/06)dé = 0 (8) 
where 0P/0M is calculated from isentropic flow re- 
lations. 

Next, if p = P/P, where P, is the total pressure for 
the flow along the jet boundary, then Eq. (8) becomes 


(0p/OM)dM + (dp/d0)do = 0 (9) 


P,is assumed constant along the boundary, since the 
waves are of infinitesimal strength. 
For isentropic flow, p is related to .V as follows: 


p= (lt (vy —-— DM/2)0-"7 (10) 
so that 
(0p/0M) dM = (0p/0M?*) dM? = 

—¥/2(p d M2/}1 + [(y — 1)/2]M4}) (11) 


From linearized theory, the change in pressure due 
to a small change in direction is 


AP = (y P M?/V M? — 1)d0 (12) 


Hence, in terms of a vanishingly small change in 
pressure and direction, 


dp/20 = (y p M2/V M2 — 1) (13) 


Substituting Eqs. (11) and (13) in Eq. (9), one ob- 
tains the following relation between 6, the angle meas- 
ured from the original direction of boundary flow at 
the nozzle exit, and 7, the Mach Number at any point 


along the boundary: 
(—y/2) (b/}1 + [Cy — 1)/2]M}) dM? + 
(y p M?2/V M2? — 1)d0 = 0 (14) 


Substituting 8 = VM? — 1 into Eq. (14) and re- 
arranging, one obtains 


d@ = [2/(y — 1)] [6?/(8? + 1)] X 
iy + D/(y — D1) + Bde (15) 
Eq. (15) can be integrated now, from the nozzle-exit 
conditions after expansion, to any point on the bound- 
ary. At the nozzle exit, 8 = 8, (\J = A/,) and @ = 0. 
Thus, 


a= 
V(y + 1)/(y — 1) tan [Wy = )/(y + 98) - 
tan-! 8B — iV (y + 1)/(y -— 1) Xx 


tan—! [(V(y = 1)/(y +1)8] —tan-!8,} (16) 
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Hence, with the given approximations, Eq. (16) indi- 
cates that for the condition of constant pressure to hold, 
the turning angle, 6, at any point, is simply the differ- 
ence between the Prandtl-Meyer angle associated with 
the Mach Number corresponding to the area ratio at 
that point and the Prandtl-Meyer angle associated 
with the Mach Number, .V/,. 

With @ known, one can calculate the jet boundary 
in terms of r and x, since if w is the angle made by the 
tangent to the jet boundary with the axis, then (Fig. 9) 


tan y = (dr/dx) (17) 
But y=a-0 (18) 


and if one defines a dimensionless y and x in terms of the 
radius of the exit section of the nozzle then Eq. (17) 
can be written as, 


d(x/ry) = d(r/ry)/tan (a — @) (19) 


Hence, x/ry at any point (where x = 0 at the nozzle 
exit) is as follows: 


r/*N 
(x/ry) = f [d(r/ry)/tan (a — @)] (20) 
1 


The integral in Eq. (20) may be calculated numerically, 
using Simpson’s rule of integration. At any r/ry, 
the corresponding area ratio (A/A*) may be calcu- 
lated, using as(A /A*), ;,,~; the value which corresponds 
to V,. Thus, 


(A/A*),+y = (r/rv)? (A/A*)u om, 


Knowing the area ratio, the Mach Number and @ can 
be found from compressible flow tables. For any 
given set of nozzle parameters, a can be calculated using 
Eq. (6). 

In Figs. 10(a), (b), and (c), plots of 7 vs. x are given 
for nozzle angles of 6y = O for various Py/P., when 
My = 1, 2, and 3. On the same graphs are plotted 
the characteristic solutions reported by Love and 
Grigsby. 

It is clear that this approximation is good only when 
the angle which the jet boundary makes with the axis 
at the nozzle lip is relatively small—i.e., at relatively low 
pressure ratios, Py/P.,. Furthermore a nozzle Mach 
Number effect exists. At Mach Number, Vy = 1.0, 
for example, the maximum percentage difference, based 
on the characteristic solution, goes from 5 per cent 
at Py/P,.. = 2, where P:/P.. = 3.8, to 9 per cent at 
Py/P. = 20, when P,;/P, 2 80. At My = 3, how- 
ever, this percentage difference is 16 per cent at Py /P. 
= 2, where P,/P.. = 74, and 29 per cent even at Py /P. 
= 10, where P,/P.. = 370. Hence, the approximation 
is limited by pressure ratio and nozzle Mach Number, 
depending on the accuracy desired. 


evar uM 


Fic. 11 (right). Comparison of approximate jet boundaries with 
the boundaries seen on a shadowgraph. Top to bottom: 
(a) My = 1, P:/Po = 18.5, Py/Po = 9.77 


(b) My = 2, 6y = 10°, P:/Po = 24.6, Py/Po = 3.14 
(c) My = 2, 6y = 10°, P;/Po = 62.2, Py/Po = 7.95 
(d) My = 2, 6y = 10°, P:/Po = 138.7, Px/Po = 17.7 
(e) My = 3, dy = 10°, P,/Po = 139.5, Py/Po = 3.80 
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To give some idea of the comparison between these 
calculations and experiment, the calculated curves 
have been drawn on the experimental photographs in 
Figs. 11(a) through (e). In each case the pressure 
ratios of both the experiment and the theoretical calcu- 
lations agreed within a few per cent. Again, at rela- 
tively low pressure ratios, the agreement is quite good, 
while at very high pressure ratios it is quite inaccurate. 

It should be noted that, with the approximate method 
of calculating the inviscid jet boundary, only the diverg- 
ing part of the flow can be considered. Thus, asy — 0, 
i.e., as the jet boundary tends to become parallel to 
the axis, x > ©, since 1/tan ~ ©. This restriction 
was pointed out also by Love’ in a recent note in which 
he used a similar method to calculate the boundary be- 
tween two supersonic streams with different Mach 
Numbers and pressures. 


CONCLUSIONS 


A method for calculating the position of the first 
normal shock in the jet flow behind highly underex- 
panded nozzles has been found. It gives good results 
when compared to experimental results, for a range of 
pressure ratios of 5 < P,/P., < 140 and for range of 


nozzle exit Mach Numbers of 1 < My < 3. 
The approximate method presented for calculating 
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the boundary of a supersonic jet exhausting into an 
atmosphere at rest is much easier to apply than a char- 
acteristic solution. However, its validity is limited by 
the accuracy desired to relatively low nozzle-exit-to- | 
atmospheric-pressure ratios, and relatively low super- | 
sonic Mach Numbers at the nozzle exit. 
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An Approximate Analysis of Nonlinear 
Flutter Problems 


S. F. SHEN* 
University of Maryland 


SUMMARY 


The approximate analysis of nonlinear flutter problems by the 
method of Kryloff and Bogoliuboff is outlined. The method is 
presented in a simple fashion based upon the concept of ‘‘har- 
monic balance.’’ Because of the nature of the response of flutter 
systems to sinusoidal excitations of various frequencies, it ap- 
pears that the method should be particularly suitable in the 
present application. A nonlinear system is reduced to an equiva- 
lent linear system, of which some of the parameters become func- 
The flutter behavior of non- 
linear systems can thus be explained, often very simply, in terms 
of the knowledge of the linear system. 

A few examples of wing-control surface flutter with nonlinear 
structural stiffnesses, solved by an analog computer by Woolston, 
The analytical results predict, at least 


tions of the amplitudes of motion. 


et al.,).? are re-examined. 
qualitatively, the behavior as described there. 

It is difficult to compare the flutter boundaries of the present 
analysis and the published analog results because of the use of 
different amplitude parameters. The agreement is good in one 
case where comparison is possible. 

The method is next applied to the finite amplitude flutter of 
buckled panels in a supersonic stream, first studied by Fung‘ 
with a two-mode representation, to illustrate its capability in 
handling cases with many degrees of freedom. Without using 
the mode approach, it is shown that the nonlinear behavior may 
be deduced from a knowledge of the linear system, which in this 
particular case becomes identical to a problem treated by Hedge- 


peth.® 
INTRODUCTION 


I A BROAD SENSE, flutter is but one type of self- 
excited oscillation complicated by the aerodynamic 
The usual linear 


21,—) theory assumes both the structural and the aerodynamic 
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properties to be independent of the amplitudes of os- 
cillation. Solution of the linear flutter problem has 
long since become a routine, if often tedious, matter. 
However, to every engineer it is probably obvious that 
icertain structural nonlinearities must always exist 

for instance, in the form of backlash between elements. 


» Aerodynamically, too, nonlinear variation with the am- 


plitudes of motion is expected to be appreciable in 
such problems as stall flutter, or as the free-stream 
Mach Number is increased. Recently, a number of 
typical flutter problems with certain nonlinear struc- 
tural properties have been investigated on an analog 
computer.’ *? The results are illuminating, but one 
cannot get a feeling of the role of the nonlinearities by 
examining a set of specialized data only. The purpose 
of this paper is to point out the convenience of an ap- 
) proximate analysis in dealing with flutter problems with 
fany kind of nonlinearity. The method is that of 
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Kryloff and Bogoliuboff,* which is well known in the 
theory of nonlinear oscillations, and has been applied 
at least once in a nonlinear flutter problem.‘ The 
simplicity of its basic concept and the flexibility in its 
applications seems to be worthy of the attention of 
aeroelasticians. 

We shall first briefly summarize the procedure of 
analysis and discuss its application to flutter problems. 
It is observed that, for problems of the type treated 
by Woolston et al.,'* the answer can be directly ob- 
tained if the behavior of the linear system with various 
stiffnesses and structural damping coefficients is known. 
Our results for some of the examples may be inter- 
preted to represent, without contradiction, the phe- 
nomena of “mild flutter,” “violent flutter,” and ‘‘di 
vergent flutter,’’ as obtained from analog computa- 
tions. Direct comparison of the flutter boundaries, 
however, is difficult because of the different choice oi 
the amplitude parameter. The amplitude of steady 
oscillation in one of the degrees of freedom is used as 
the reference parameter in this analysis, while the initial 
amplitude of disturbance, applied to a chosen degree 
of freedom, was adopted in references | and 2. In one 
example, the limiting amplitude was given in the analog 
record, and the agreement of the flutter boundary with 
ours is quite good.'*. 1* The technique is next applied 
to the problem of flutter of a buckled plate in a super- 
sonic stream, in a way different from Fung’s original 
approach.‘ By avoiding the use of a limited number 
of modes, our treatment should be more satisfactory. 
Although the numerical work does not seem to be much 
more tedious, it has not been carried out in this paper. 
We only note that the problem is then reduced to one 
mathematically identical to the linear case studied by 
Hedgepeth.® 


THE METHOD OF KRYLOFF AND BOGOLIUBOFF 


We shall present only a simplified version of the 
method of Kryloff and Bogoliuboff. Furthermore, we 
restrict ourselves to their lowest approximation and 
assume a heuristic point of view. Readers interested 
in the general theory and more rigorous discussions are 
referred to the original work (e.g., see reference 3). 

We begin by an example. Consider the free oscilla- 
tion of a one degree of freedom case with a cubic spring, 


#+ F(x) =0 ) 
F(x) ax + Bxf 


Let us assume that a periodic 


(1) 


a and 8 being constants. 
solution exists, in the form 
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Fic. 1. Flutter coefficient Vr/wabr against frequency ratio 
w8/wa With and without damping. 





> a, sin nwt (2) 


n=1 
Substituting Eq. (2) into Eq. (1), we get 
>> (—n’w’a, + b,) sin not = 0 
2=1 


F(x) ~ » a b, sin nwt 


n=1 } 


with 
To balance the fundamental harmonic (n = 1), there 
must be 


—aw* + db; = 0 | 


i.e., w? = b,/a; = 2/ray) f F(x) sin wtd(wt) ( (4) 
0 


If the fundamental harmonic dominates, then 
T 
w? > (2/ma;) f F(a, sin wt) sin wld (wt) (5) 
0 


or 
w Sat (3/4)Ba;? (6) 


This is the same answer as obtained by other methods 
for the present problem (‘‘Duffing Equation,” e.g., 
see Stoker,® p. 85). The result is valid as long as the 
assumed property holds true—namely, that the system 
admits a periodic solution dominated by the fundamen- 
tal harmonic. The general application is then to as- 
sume a periodic solution and balance only the funda- 
mental harmonic. For this reason, the approach has 
been called by Kryloff and Bogoliuboff the method of 
“harmonic linearization”’ or ‘‘harmonic balance.” 
Although originally devised for one degree of freedom 
systems (second-order differential equations) with a 
small parameter for the nonlinearity, the method ap- 
parently has been developed and widely applied in 
problems of automatic control with ‘strong’ nonlin- 
earities in the U.S.S.R. To justify the application, a 


recent paper by Popov’ is concerned with the quite 
general differential equation: 


AERO/SPACE 


SCIENCES—JANUARY, 1959 


O(p)x + R(p) F(x, px) = Ol 


~] 
ae 


p=d/dt f 


where Q(p) and R(p) are polynomials and F(x, px) a 
nonlinear function. Assuming again that an approxi- 
mately sinusoidal solution exists, 


a; sin wt + ey(t), e<] 


e 
I 


> a, sin (nwt + ¢,) 


n=2 


< 
I 


one can proceed as in the simple example above by ex- 
panding F(x, px) also into Fourier components, 


F(x, px) ~ F(a, sin wt, ayw cos wt) + 
€ >. c, sin (nwt + 0,) ; 
n=0 (9) 
F(a; sin wt, dy cos wt) ~ >. b, sin (nwt + yp) 
n=0 


By substituting Eqs. (8) and (9) into Eq. (7), therej 


follows 
Q(p) a sin wt + 
eO(p)y + R(p) zu b, sin (nwt + Wp) 





eR(p) >> Cn sin (nwt + 0,) = 0 (10) 


n=0 


a REG 


If we neglect O(e), and balance only the fundamental} 
harmonic, 


aQ(p) sin wt + b,R(p) sin (at + yi) = 0 (11) 


Since b; and y; both depend on a, Eq. (11) leads to two! 
relations which may be solved for w(a;) and y(a;).* 
Thus, 


a => b,| R(iw)/Q(iw) | ( 19) 
wSr- arg { R(iw)/Q(iw)} § ( 7 
Proceeding to balance the higher harmonics, Popov) 
further arrived at the following condition on the . system 
and the nonlinear function: ; 
DPE » R(inw)? <b; R(ia)|? aa} 
> Q(inw) Q(tw) 
in order that the procedure may be justified—i.e., 
¢ < 1—in the presence of ‘‘strong’’ nonlinearities. 
Physically, one can think of the condition as requiring) 
the system to possess a “‘filtering’’ property, which) 
filters out most of the higher harmonics that may come) 
out as the output of the nonlinear element before the’ 
signal returns as the input to it again. 
It should be mentioned that rather similar ideas have) 
been aatepemtontty developed in the United States) 


SNE ae 


. if R(p) has a constant term, then do) sin Yo must be reroll 
otherwise a steady-state deviation will follow, contradicting they 
assumed periodic motion. In terms of the property of the non! 
linear function F, there must be, therefore, 


2 
f, F(a sin wt, dyw cos wt)d(wt) = 0 
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APPROXIMATE ANALYSIS OF 

although at much later dates, also in the field of control 
engineering. The term ‘‘method of describing func- 
tions’’ was offered by Kochenburger,’ for a procedure 
which is essentially identical to the above. An im- 
provement for the next higher approximation has been 
suggested by Johnson,’ who used the term “‘sinusoidal 


analysis.” 


APPLICATION TO FLUTTER WITH STRUCTURAL 
NONLINEARITIES 


Let us now apply the method of ‘harmonic balance”’ 
described above to the cases studied by Woolston, 
et al.!) * on an analog computer. The bending-torsion- 
control surface flutter of a two-dimensional airfoil-flap 
system in incompressible flow, with linear structural 
restoring forces, is described by the following equations 
in standard notations: 


mh + Sai + S83 + ma,2h = —L (14) 
Sh + Ia + [(c — a)Sg+ I,]8 + /,w,’a = M, (15) 
Ssh + [(c —_ a)Sz + Ig \& + 1,8 + T3378 = I (16) 


The simple types of structural nonlinearities treated in 
references | and 2 all amount to replacing the stiffness 
terms in Eqs. (14), (15), and (16) by nonlinear functions 
K, (hk), K,(a), and A,(@), respectively. 

Suppose now we look for periodic solutions of the 


nonlinear cases, of the type 


t 
h ~ a, at aye, B ~ Boe 


According to our method, we feed these forms into A,, 
K,, and Kg and retain only the fundamental harmonics 


which come out of the latter. Thus, 


K. ~ Tq’ 200K (a0), 


Ke — T gue’ *BoK g’ (Bo) 


K, ~ mo,/*hoK,,' (ho), 
(17) 


where w,’, w,’, and wg’ are reference values for the re- 
spective frequencies and K,’, K,’, and K,’ represent the 
effects of nonlinear stiffnesses depending on the respec- 
tive displacements. One may then solve the flutter 
determinant in the usual that 
oscillation amplitudes (/, ao, 89) appear as parameters, 
thus providing the effects of nonlinearity. 

For abbreviation, we shall call K,’, K,’, Ag’ the 
“describing functions’ after Kochenburger.* When 
they are real, we have simply a set of different, equiva- 


manner, except 


lent linear frequencies w,”, w,”, ws”, depending upon the 
amplitudes. By specifying the values of ho, ao, Bo, one 
can immediately look up the flutter speed and 
frequency from survey plots such as provided by Theo- 
dorsen and Garrick.” Fig. 1 is reproduced from their 
report as an illustration. The linear stiffness (for in- 
finitesimal amplitudes) is assumed to be at point A, 
the nonlinear stiffness (for a set of specified amplitudes) 
is assumed to be at point B. The changes of flutter 
speed and frequency from those corresponding to a 
linear system are perfectly clear. 

If the ‘‘describing functions,” K,’, K,’, and Kg’, are 
complex, our equivalent linear systems will have struc- 
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Fic. 2. (a) top. Backlash nonlinearity. (b) bottom. Charac- 


teristics of spring with free play. 


tural damping coefficients which also depend on the 
amplitudes. One must then look for the nonlinear 
flutter speed and frequency from the survey plots which 
include the structural damping coefficients (g,, g., 
£3) as parameters. There are only a few plots in refer- 
ence 10 which show the effect of the structural damping, 
but to provide for these plots is a routine matter. 

To be more concrete, we use the backlash as an illus- 




















tration. Suppose the control surface restoring moment 
versus displacement is as sketched in Fig. 2(a). If the 
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Fic. 3. (a)top. Hysteresis nonlinearity. (b) bottom. Charac- 


teristics of spring with hysteresis. 
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Fic. 4. (a) top. Backlash, preload included. (b) bottom. 
Characteristics of spring with free play, preload included. 


displacement {8 is taken to be real, clearly A, is also 


real. The output (moment) due to a sinusoidal input 
is shown in Fig. 2(b). Thus, 

K, = 0, 0<at<¢ 

Kz, = K£(sin wo — sing), ¢<wt<r-—¢ (18) 

Kg = 0, r—G<wot<r+¢ 
etc., with @ = sin—!(6/28)) < w/2, K = constant. 
Hence, for the fundamental harmonic, 

Kg = KBo[1 — (26/m) — (sin 26/7) | (19) 

As the amplitude 6» increases, ¢ > 0, Kg > Af. On 


the other hand, if 8) < 6/2, then Ag = 0; 1.e., the con- 
trol surface is free of restraints. The variation of flutter 
speed and frequency is easily traced out between these 
two extremes as a function of 0/6. 

An example of complex A, is provided by the cases 
of hysteresis, a special case being shown in Fig. 3. 
Taking fy) real, we have the output of the nonlinear ele- 


ment, 


K, = — (H/2) + (1/i)hocosat, O< at < 1/2 | 
K, = —(H/2), 7/2 < ot < (7/2) +¢ 
K, = (H/2) — (1/5)ho cos wt, (2/2) + 


o < wt < (3/2)r (20) 


K, = H/2, (3/2)r < ot < (3/2)r + @ 
K, = — (11/2) + (H/é)ho cos wt, (3/2) + 
o< wt <2n | 


The first harmonic now 


with @ = sin—!(6/ho) < 7/2. 
consists of both sine and cosine components: 


SCIENCES—JANUARY, 1959 
K, ~ {—(2H cos $/m) + (H/6)ho[1 — (6/r) + 
(sin 26/m)]} cos wt + [—(2H sin ¢/r) + 
(H/6)ho(sin 26/7)] sin wt (21) 
or, in complex notation, 


K, — Kei —y) \ 
> (22) 
, r 9 x 9\1/2 , , / 

K = (K+ &,*)"", ¥ = tan“(K, K.)§ 
where K, and K, are the coefficients of cos wf and sin 
wt, respectively, in Eq. (21). The correspondence to a 

structural damping coefficient is obvious. 

A further example of interest is the nonlinear spring 


} 
ii 
; 





which offers stiffer resistance to movements in one di- | 


rection than the other. 
for a spring with backlash and preload, the latter shift- 
ing the flat part of the diagram to the right of the origin. 
It is easy to see that a constant term will be generated 
by the nonlinear spring when a pure sinusoid is im- 
pressed on it. 
26 is not satisfied. We must, therefore, allow a steady 
deviation in addition to the fundamental harmonic: 


a= A+ asin of 
Then the output is 


K, = K(A + asinwt), x<ow<¢, 
wrt+tx—-d@<wt<2r +x 

K, = K(A + asin gd), ¢< wt < y, 
e™+tx-v<at<r+x-¢ 
K(A + ao sin wt — ap sin Y + apsin >), 
Y<w<r-—-ytx 


K, = 
with @ = sin—'(6p/ao) < 7/2, 
y = sin=! [(6p + 67)/ao] < 


sin~!(—A/apo). J 


Jé 
w/2, 
x = 


The Fourier analysis of Eq. (23) follows immediately. 
The steady-state deviation A is to be determined in the 
complete equation of motion by balancing the zeroth 
harmonic. In this particular example, we note that 
the fundamental harmonic happens to be identical 
with what would have been the result without allowing 
for the steady deviation A. 

For fundamental harmonic only, the nonlinear spring 
with backlash and preload thus is represented by 


| fiat K, ay < dp, ) 
K, ~ K{((1/2) + (¢/r) + (sin 2¢/27)], 
Op <a< bp + 67, 


, ha o sin 2¢ f 

K,~ K E +-+ ie | 
“ T Tv Tv 

2 cos y [(3/2) sin y — sin 9] 


T 


| ay > Op + Or | 


We do not attempt here to justify our application to 
flutter problems by checking with Popov’s condition 
[Eq. (13)]. The quantity Q(7w) in flutter problems will 
not be a polynomial because of the transcendental Theo- 
dorsen’s function. Physically more important, however, 


In Fig. 4 we show the diagram | 


Hence, the condition in the footnote of p. | 





(23) 
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is the fact that flutter systems generally do have a very 
strong “‘filtering’’ property. Let us assume that we 
have determined the equivalent linear system by balanc- 
ing only the fundamental harmonics as described above. 
We may now replace the actual nonlinear system by the 
equivalent linear system p/us an element which gener- 
ates the higher harmonics. As sketched in Fig. 5, let 
the input be the fundamental harmonic of the flutter 
frequency. Then it reaches the input end of the non- 
linear element with the same frequency, after going 
through the equivalent linear system. Now, of course, 
higher harmonics are generated, which will again come 
back to the input end of the nonlinear element as the 
response of the equivalent linear system to the higher 
frequency input. It is a general characteristic of linear 
flutter systems that the response to sinusoidal input is 
sharply peaked only at the flutter frequency. We 
are, therefore, in little danger that any appreciable 
amount of higher harmonics may persist even for our 
nonlinear problem. This argument suggests that the 
method should be particularly suitable in flutter prob- 


lems. 


NUMERICAL EXAMPLES 


We shall now apply the above procedure to analyze 
a few examples which are covered in references | and 2. 
The properties of the systems are summarized in Table 
l. 
Case (1)—-Flutter of a Two Degree of Freedom System 

With Free Play 

Both experimental and analog results are given in 
The degrees of freedom are pitching and 
The system is linear in translation but 


reference |. 
translation. 
has a nonlinear torsional stiffness as shown in Fig. 4. 
In reference | it is mentioned that a certain preload 
(i.e., dp in Fig. 4) was incorporated in both the analog 
and the experiment, but the exact amount is not given. 
It turns out to be variable (between 0.01° and 0.08°), 
depending upon the ratio V/Viin. (see reference 13). 

In our calculation a curve of flutter velocity versus 
K,,’ is first prepared (see Fig. 6). The describing func- 
tion K,’ is obtained from Eq. (2+) as 




















a, Go 


Hl} 1 


EL: Equivalent linear system 

N : The nonlinear element 

q, : Output of N, input of EL 
G2: Output of EL, input of N 


Equivalent linear system. 
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Fic. 6. Flutter coefficient Vr/wab against equivalent spring 
constant, 
x” = |]. co < bp 
K,’ = (1/2) + (¢/r) + (sin 26/27), 


Op < ay < Op + br 
K,’ = (1/2) + (6/7) + (sin 26/7) — 
») 


(W/m) —} } 


cos ¥[(3/2) sin y — sin o|/7j 
ay > Op + Or 


l(6p/ao) & w/2, 


y = sin 


with @ = sin 
(6p + dr) ay | . w/2 

The reierence frequency w,’ is taken to be the linear 
value (6¢ = 0, K,’ = 1). ’ with 
a tor fixed 6p and 6r7 is obvious from Fig. 4 when ay > 


The variation of AK, 


6); it first decreases toward a minimum, but eventually 
rises again toward 1. 

It is then a simple matter to pick out the flutter 
speed variation versus ay by combining Eq. (25) and 
The results for 67 = 0.5° and two values of 
The region to 


Fig. 6. 
dp (0° and 0.2°) are plotted in Fig. 7. 
the left of each curve is stable; to the right, unstable. 
There is a minimum speed, according to our theory, be- 
low which the system is stable with regard to all ampli- 
tudes ay. Immediately above the minimum speed, 
stable oscillations of amplitudes of about 1° to 2° would 
This should correspond to the so-called ‘‘mild 
For higher speeds, the 


occur. 
flutter’ in references 1 and 2. 
oscillations should still be stable but with rapidly in- 


TABLE | 
Properties of Systems 


Three Degrees 
of Freedom 


Two I Jegrees 


Parameter of Freedom 


a —(). 406 —().2 
Xa 0.042 0 

ra? 0.1603 0.5 
w, (rad. /sec. ) 58.94 81.7 
wa (rad./sec.) 81.24 136.0 
wp (rad./sec. ) Variable 
€ 0.6 
xB 0 

rp? 0.002 
1/k 156.405 12.0 

b (ft.) 0.5 2.375 
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Bending-torsion flutter with free play and preload in 
the torsional stiffness. 


Fic. 7. 


creasing amplitudes, corresponding to the “violent 
flutter’’ of references | and 2. The analog flutter 
boundary does not resemble our result, but there the 
amplitude of the initial disturbance which excited the 
flutter motion was used as the reference parameter in 
place of our steady amplitude at flutter. 

While direct comparison is difficult, our curves are 
believed to be at least qualitatively correct. In terms 
of the effective stiffness when the amplitude covers at 
least part of the flat spot, it should be clear that the 
weakest torsional stiffness occurs at ay = 6p + 5dr. 
Thus, the calculated minimum flutter speeds in Fig. 7 
are at ay = 0.5° and 0.7°, respectively, for the assumed 
6p’s. For much larger ao’s, the behavior approaches 
the linear curve as it should.* 


Case (2)—Wing-Control Surface Flutter With 

Free Play and Preload 

This is a three degree of freedom problem, the proper- 
ties of the system being listed in Table 1. We assume, 
as in reference 2, that the restoring moment to the con- 
trol surface deflection is represented again by the dia- 
At small amplitudes the system is 
In fact, 


gram of Fig. 4. 
linear and has been analyzed in reference 9. 
we shall directly use two survey plots for g, = 0, 
£4 = 0.02, reproduced as Fig. 1. 

Except for replacing a by 8, the nonlinear element is 
again as described by Eq. (25). Therefore, we easily 
obtain the curves in Fig. 8 by a combined use of Eq. 


* The flutter boundaries of the present theory were inadvert- 
ently compared against those of references 1 and 2, by Shen and 
Hsu in reference 12, without their noticing the different choice 
of the amplitude parameter. Objection was raised by Woolston 
and Andrews (see reference 13), who made available in the same 
reference some information of their values of 5p (varying between 
0.01° and 0.08°, depending upon V/Vjin.) and ao for the case in 
Fig. 7. The dotted curve constructed on the basis of these data 
is the only one which could be compared with our theoretical 
curves (see reference 14). The agreement shown in Fig. 7 is 
quite reasonable 
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i 
: 


1, assuming two values 6p/éir = 1 and 
According to our analysis, for small 


(25) and Fig. 
bp/br = 0:5. 
values of 8)/67, the system is unstable for V/Viin. >| 
1.0, but a limit amplitude 8/67 will be maintained if the 
ratio V/Viin. is not much greater than 1.0. But at] 
these same speeds, if the amplitude is increased, a} 
boundary for divergent flutter will be reached, and, | 
above a maximum flutter velocity, divergent flutter} 


~ > emeevermte 


always occurs. 

The qualitative behavior agrees with the analog re- 
sult reported in reference 2. The flutter boundary 
sketched in reference 2 is of a different nature, but, 
again, the amplitude parameter used was different.*| 
BUCKLED PANEL IN A SUPERSONIC 
STREAM 


FLUTTER OF A 


To illustrate the flexibility of application of the] 
method, we next consider the nonlinear flutter problem 
of a two-dimensional buckled panel in 
stream, a system with an infinite number of degrees of 
freedom. The problem of the buckled plate having an 
axial load proportional to the change of length of the 


supersonic 


plate, and assuming quasi-steady aerodynamic force, 
was first studied by Fung.‘ In Fung’s notation, the 


following nonlinear equation results: 


(O*Y/OX*) + (1 + x)(0?Y/0X?) = 
—Q[x(OV/OX) + a(OV/dr) + ¥(02¥/0r2)] 


(26) 
°n 
x = A? — (2/r) | (OY /OX )*dX 
The boundary conditions for a simply supported panel 
are 


"sab OX?) x <0 == 
d°V/0X2)x-, = 0 (27)| 


Y(0, r) = Y(a, 7) = 0, 


We shall only mention that ),°, Q, a, y are all constant 
parameters, and that the flutter problem is reduced 
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in the aileron stiffness. 
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APPROXIMATE ANALYSIS OF 


mathematically to the determination of the eigenvalue 
Q, for assigned values of ),*, a, and y, when neutral 
oscillations become possible. The nonlinear term in 
Eq. (26) is x(0?Y/0X?), causing the solution to depend 
upon the amplitude of the oscillation. The reader is 
referred to Fung’s original work for the physical mean- 
ings of these symbols. 

Fung’s analysis was carried out on the basis of two 


“modes.” By putting 


Y ~ B,(r) sin X + B.(r) sin 2X 


Static equilibrium was 
The dynamic 


the Galerkin method was used. 
found by taking B, and B, as constants. 
stability of these static equilibrium configurations were 
discussed first on the basis of small perturbations, then 
by adapting the Kryloff-Bogoliuboff theory. It may 
be noted that the results for the finite disturbances were 
represented in a group of charts, which are not easy to 
interpret if only the physical system is given. 

It is felt that perhaps the picture would be in better 
focus if the mode approach is abandoned right from the 
A possible source of error—that of an in- 
will also be eliminated. 


beginning. 
sufficient number of modes 
The role of flutter amplitude alone may thus be better 
assessed, as compared to the role of initial warping, for 
instance, which was treated by Fung in a later paper."! 


Let us start with the static equilibrium. Putting 
Y = f(X) 
we get, from Eq. (26), 
fo" ++ xf’ + Onf’ = °| 
(28) 


x=rAZ2—- (2 nf frdX 
0 


with boundary conditions, 


f(0) = f(r) = 0, f’(0) =f'(r) = 0 (29) 


It is obvious that, for any f(X), the parameter x is a 
constant only, and Eq. (28) remains an ordinary difter- 
ential equation with constant coefficients. By assum- 
ing the constant Q, the eigenvalue x and the eigen func- 
tion ¢(X) may be determined in the usual way. If our 
solution is to be f = Ad¢(X), the amplitude A may be 
finally determined: 


/ 


x 1/2 
A= | - x) /(2 af ordX | (30) 
0 


The functional dependence of Qcr, (the critical value of 
Q) on A is thus obtained in an inverse manner. 

Now, for the dynamic problem, we apply the concept 
of “harmonic balancing’’ by assuming again a simple 


harmonic motion. Putting 


Y = RI [g(X)e"*"] 


we get 
x =dAZ2 — (2 * fe {RI [g’(X)e*"]} dX 
” » (31) 
= 2? — (2/m)[J; cos? wr — 2]. sin wr X | 


cos wr — J; sin? wr] 
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where 


I -{ [gn (X)]*dX 
0 


Is -{ &r' (X)g,'(X)dX ? (: 
0 


I; -{ [gr’(X) }*dX 
0 


and gr, g, being the real and imaginary parts of g(X). 
As before, the term (1 + x)(0?V/O0X*) may be expanded 
into a Fourier series for the assumed Y, and only the 


Ow 
to 


fundamental harmonic will be retained: 


(1 + x)(0?Y/0X2) ~ RI [F,(X)e"] _ ) 
F, = ag’’(X) + bg’’(X) 


a = —(1/2m)[2/, + + (33) 
(1 + a)(2 + Is) + 1+ 12] | 


and g is the conjugate of g. Hence, Eq. (26) is reduced 


to 


"oe + ag” — bg”’ a Org’ - Ag = O) 
A = yQw?[1 — (ta/yw) | 


The boundary conditions remain, 
2g(0) = g’’"(0) = 0, g’’(0) = g(r) = O 


It is seen that, as A is complex for a ¥ 0, the eigen- 
function g is generally complex. 

However, for a = 0 (which occurs at Mach Number 
1/2) or |a/yw| < 1,* so that the a term is negligible, 
both A and g may be real. With real g (‘‘standing 
wave’), the coefficients a and } are real, and we may 
try to find the real eigenvalue A. Should real A be- 
come impossible, it is easily concluded that, because 
both A and its conjugate are eigenvalues, divergent 


motion must result. The equation is now 


g’’"’ + (a + b)g’”’ + Que’ — Ag = 0| (35 
P oo) 
A => 7Qw? f . 


Again (a + 5) isa constant, although depending on the 
By assuming a value for (a + 6)—i.e., an 
Eq. (35) becomes mathematically 


amplitude. 
average axial load 
identical to the problem studied by Hedgepeth.® At 
Q = 0, A gives the vibration frequencies im vacuo. As 
Q increases, a critical value Q.,. will be reached, beyond 
which at least some A’s become complex, and the mode 
o(X, Qer.) is thus determined. Both Q.,. and @(X, 
Qer.) depend upon (a + b). Let g = Ag; then (a + 
b) is a function of the amplitude A and the mode ¢ 
Weare, therefore, again in a position to define implicitly 
Qcr. versus A, the flutter boundary versus the amplitude 
of motion. It is straightforward to use Hedgepeth’s 
eigenvalues and eigenfunctions for his equation, which 


* The condition |a/yw | < 1 is often satisfied. In Fung’s 
special case, this was the basis for Eq. (47), reference 11. At 
high Mach Numbers, its value varies as M~'; see, e.g., Hedge- 


peth.® 
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Fic. 9. Stability boundaries for a plate which is subjected to 


compressive load greater than the first buckling load, yet remain- 
ing flat, in a supersonic stream 


is identical to Eq. (35), to illustrate our procedure as 
described above.f 

Without much calculation, we may draw some quali- 
tative conclusions by borrowing Hedgepeth’s data to 
illustrate the dynamic stability around a static equilib- 


rium configuration. Putting 


Y = f(X) + RI[y(X)e“"] 


we obtain, from Eqs. (26) and (28) for small perturba- 
tions, 


en 


yl’”’ + (1 of xo) yp’ 4+ (4 w)f’’ | Sf ydX 4 


0 


Gry’ — Av =0 } (36) 
jo = A, = (2 ") | f?dxX 


0 ) 
with ¥(0) = ¥(r) = 0, W’’(0) = W’'(r) = 0. When 
Eq. (36) has a solution, the panel may depart from its 
static equilibrium configuration and execute an oscilla- 
tion. But for Q > Q.,. some complex A’s occur, and the 
motion becomes dynamically unstable. 

A particularly simple case is that of f = 0; 1.e., the 


panel being blown flat in static equilibrium. Under 
this assumption, Eq. (36) is reduced to 
vr" + (LE MY” + Ory! — AV = OL 
of 


¥(0) = V(r) =0, W’(0) = v(x”) = OF 


Obviously Eq. (37) is mathematically identical to Eq. 

(34), so that Hedgepeth’s numerical data may be used. 

After proper scaling, we obtain the following: 
Mi 2 V/3 


1 ( 
Qer. 0 0.60 1.95 


— 


bo 


72 


In Fig. 9, we have replotted Fung’s Q.,. versus \; for 
g I g 


{ Calculations in this direction have now been made. The 
results will be reported separately. 


JANUARY, 1959 


dynamic stability (R = 0 curve of Fig. 2, reference 4), 
based on a two-mode approximation, and the above 
numbers deduced from Hedgepeth’s data are included 
for comparison. Our ‘exact’? solution appears to 
give 20 per cent higher Q.,. for A; < 1. 

Qualitatively, we now go back to discuss the nonlinear 
problem of Eq. (35). Because of the definitions in 
Eqs. (32) and (33) for real g, 


a+tb=14+A/ -— (3 2m) f g’*dX 
0 
=1+A),° -— (3 2n)A? | o’*dX 
0 


The effect of amplitude A is thus to correspond to a 
lower );” of Eq. (37), consequently increasing Q.,, One 
may predict, therefore, that for finite amplitude oscil- 
lations with a given initial \,”, the flutter speed ought 
to be higher than that predicted by Eq. (37). 


CONCLUSIONS 

In this paper, the analysis of nonlinear flutter prob- 
lems by the method of Kryloff and Bogoliuboff is sug- 
gested. The method is presented in a simplified form, 
according to the concept of “harmonic balancing,’’ so 
that it should demand no undue effort on any practicing 
aeroelastician. Because oi the very sharp peaking of 
the response of a system near the flutter conditions to 
forcing functions of the proper frequency only, it ap- 
pears that the method should be particularly successful 
in the present application. 

The cases of nonlinear structural stiffnesses of any 
type are the simplest to treat. According to the 
method, we only need to consider an equivalent linear 
system whose stiffnesses depend upon the amplitudes 
of displacement. When the nonlinearity is of the hys- 
teresis type, the equivalent structural damping coef- 
ficient is likewise a function of the amplitude. By 
such equivalence, the problems of references 1 and 2, 
solved on an analog computer, are handled with great 
ease. Our analysis yields results which agree qualita- 
tively with the behavior reported by references | and 
2. A closer comparison is difficult because the analog 
data were presented by using a different amplitude 
meter. In one case the comparison was possible and 
the agreement good (see reference 14). 

One advantage of the analytical method described 
in this paper is its flexibility in covering various non- 
linearities. The simultaneous existence of nonlinear 
structural stiffnesses in all degrees of freedom offers no 
more difficulty of analysis. To illustrate, we have ap- 
plied the method to Fung’s problem of a buckled panel 
in a supersonic stream in a way different from his orig- 
inal approach, which was based upon only two modes. 
Our approach admits an infinite number of degrees of 
freedom, and the numerical part does not seem to be 
much more tedious. It goes without saying that the 
very same method is directly applicable when non- 
linear aerodynamic effects might become important 
for instance, as higher Mach Numbers are encoun- 


(Continued on page 45) 
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A Comparison of Theoretical and 
Experimental Loads on the B-47 
Resulting From Discrete Vertical Gusts 


CHARLES E. JACKSON* ann JOHN E. WHERRY** 
Boeing Airplane Company 


SUMMARY 


During the past several years there have been published a 
number of excellent papers on dynamic gust analyses. These 
papers have dealt largely with the qualitative aspects of the 
subject in order to point out the significance of the dynamic loads 
as well as to determine important factors affecting them. How- 
ever, in the light of current specifications that will require dy- 
namic gust analyses of a large percentage of future designs, there 
exists a need for comprehensive quantitative comparisons be- 
tween theoretical and experimental data in order that the dy- 
namics engineer may assess current methods and future analysis 
needs. This paper presents such a study made for discrete verti- 
cal gusts on the B-47. 

The study undertaken to determine 
realistic gust loads could be predicted using presently available 
dynamic analysis techniques. Experimental gust data were ob- 
tained from a B-47 (instrumented to record the time histories of 
shears and moments at various wing stations, as well as the time 
history of the c.g. acceleration) flown over a “gust course.’ 
Theoretical dynamic analyses were made for experimental con- 


was whether or not 


ditions in which discrete gusts were encountered with a rela- 
The results of these analyses are com- 
The experimental gust 


tively smooth entry. 
pared with the experimental results. 
loads are also compared with experimental maneuvering loads in 
order to evaluate the accuracy of the gust loads formula used in 
conjunction with a steady-state analysis. 

In initial studies the effect of various analysis techniques were 
investigated. Results are presented which show the effect of 
including various airplane elastic modes and of computing loads 
by an approximate method. 

The equations of motion and the finite difference equations 
are derived and presented with the load equations in the Appen- 
The equations are in a form suitable for the solution of large 
Final results are 


dix. 
problems by digital computing techniques. 
given in the Conclusions. 


SYMBOLS 


rectangular matrix 


} } = column matrix 
| = row matrix 

[°] = diagonal matrix 

Hh = matrix of inertia load coefficients (load summation 
method ) 

[D| = matrix of air load coefficients (load summation 
method ) 

{71} = matrix of dimensionless quarter chord displacements 
in terms of the normal mode spatial data 

2 = reference inertia, /, = m,b,? 

[J] = matrix of dimensionless mass data 

K(S) = Wagner lift growth function 

|M} = matrix of vertical bending moment 
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Mach Number 
acceleration at the airplane c 


=} 


matrix of generalized forces 

dimensionless time = Vt/b, 

induction matrix 

matrix of torsional moments 

vertical gust velocity 

matrix of vertical shears 

airplane speed 

work done on the airplane by the externat air forces 

two-dimensional lift curve slope at M = 0 

local semichord 

matrix of bending moment coefficients (mode dis 
placement method ) 

reference semichord 

gust length 

gravitational acceleration 

vertical displacement of the quarter chord referred to 
the local semichord 

lift at the quarter chord on a streamwise surface 
element 

pitching moment the 
streamwise surface element 


about quarter chord on a 

reference mass 

number of normal vibration modes 

pressure 

matrix of generalized coordinates 

lag in gust entry time for the 7th air load panel due to 
wing sweep and horizontal stabilizer location 

matrix of shear coefficients (mode displacement 
method ) 

real time 

time the leading edge of the 7th airload station has 
been in the gust 

matrix of torsion coefficients (mode displacement 
method ) 

prefix, indicates an incremental change 

sweep angle 

variable of integration over an interval of real time 

summation 

matrix of airplane acceleration coefficients 

Kussner lift growth function 

angle of attack 

prefix, indicates a variation 

downwash angie at the horizontal stabilizer 

dimensionless displacement at the inertia reference 
stations 

pitch angle 

pitch rate 

distance measured aft from the c.g. of the airplane 
to the leading edge of the wing or horizontal sta- 
bilizer at the 7th air load station 

air density 

variable of integration over an interval of dimension- 
less time 

matrix of dimensionless normal mode spatial data 

natural frequency in the 7th normal vibration mode 


reference frequency = V/b, 
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INTRODUCTION 


fbn TREND toward designs which have inherently 
large structural dynamic effects due to gust en- 
counters is common knowledge. This trend has cul- 
minated in both civil and military requirements that 
dynamic gust analyses be made to determine gust 
loads and to substantiate airplane safety.'? While a 
number of papers have been published dealing with 
various aspects of the problem,*~> data which can pro- 
vide a basis for assessing the accuracy of present anal- 
ysis techniques in producing load estimates suitable 
for design have not been generally available. The B-47 
gust study provided an excellent opportunity to obtain 
such data for discrete vertical gust encounters. 

Although, as several authors have pointed out, the 
discrete gust concept does not provide an accurate 
model of atmospheric turbulence, comparative studies 
such as the present one are useful because (1) some gust 
loadings can be well represented as discrete pulses— 
blast loads, for example, are of this type—and (2) dis- 
crete gust studies provide a relatively simple method 
of evaluating analysis accuracy. 

While this study was undertaken for a specific air- 
plane, it is felt that information gained from it is gen- 
erally applicable to gust analyses of other flexible air- 
planes, and will be useful in establishing techniques for 
analyses of future airplanes. 


EXPERIMENTAL PROGRAM 


As part of a general experimental flight gust study, a 
series of records were taken at low altitudes in an at- 
tempt to obtain data from discrete vertical gust en- 


1959 
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counters. Runs at speeds from 225 knots to 400 knots 
were made at an absolute altitude of about 150 ft. over 
the crest of a saddle formation in the Wichita Moun- 
tains of Oklahoma. The pressure altitude during these 
flights was about 2,000 ft. 

The test plan called for flights at about 125,000 Ibs. 
with both forward and aft c.g.’s and runs at about 
180,000 Ibs. with a single c.g. location. The weight 
variations experienced during the high weight runs 
were so large that the gust data were relatively useless 
for the present study. The data taken at the low 
gross weights provided some excellent samples of dis- 
crete gust encounters, however, and these formed the 
basis for the present comparative study. 

The test airplane was instrumented extensively, as it 
had been used previously in the B-47 flight loads sur- 
vey and structural demonstration programs. A dia- 
gram of instrumentation significant to this paper is 
shown in Fig. 1. In addition to that shown, instru- 
mentation to indicate control surface deflections was 
included and proved to be valuable in helping to dis- 
tinguish between gust and maneuver flight conditions. 
Gust velocities were recorded by a differential pressure 
probe and an angle of attack vane. Both of these 
pickups were mounted on a boom from the aircraft 
nose. 

In analyzing the experimental data, only conditions 
for encounters preceded by what appeared to be a rela- 
tively steady condition, or data for gusts so large that 
the effect of the condition prior to encounter were 
negligible, were used. The wing loads for the gust 
data reduced showed negligible antisymmetrical com- 
ponents. This was evidently due to the fact that only 
gusts which were approximately symmetrical resulted 
in appreciable wing loads. A typical record is shown 
in Fig. 2. 

The onset of the gust was determined from the time 
at which the c.g. acceleration departed from the pre- 
gust reference. All the increments were measured 
from data values at that time. The gust length was 
estimated as the distance traveled during the time re- 
quired for the c.g. acceleration to complete the first 
half cycle of vibration. Theoretical analyses, for the 
range of gust lengths encountered, had previously shown 
that this assumption was sufficiently accurate for the 
present study. 

For the peak load studies the load data were reduced 
toa ratio, 


(A load) peak, '(AN) peak 


Generally these peaks did not coincide timewise. 
The incremental loads per unit c.g. acceleration were 
used to eliminate the complication of trying to account 
for discrepancies between actual and apparent gust 
velocities. Estimates of the gust velocities were made 
using both gust probe and gust vane data, however, 
in order to provide an added check of the theory. The 
gust velocities were estimated with the following ap- 
proximate formula: 
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THEORETICAL AND 


U, = V tan Aa + gf AN(r)dr 


The integral was taken over the interval from the gust 
encounter to its peak value. Obviously the influence 
of pitching has been neglected. Data recorded during 
the test substantiate this assumption, however. 


THEORETICAL ANALYSIS 


The theoretical analyses were the chosen mode type. 
Rigid body pitch and translation and a number of air- 
plane vibration modes were used as the generalized 
coordinates. The vibration modes were computed 
from predicted stiffness and mass distribution data 
and are listed in Table 1. 

The derivation of the equations of motion, their con- 
version to finite difference equations, and the load 
equations are presented in the Appendix. The theo- 
retical analysis technique is an extension of unpub- 
lished work by K. R. Thorson and R. B. Skoog of Boeing 
Airplane Company. 

Significant features of the analysis are: 


Structural 

(1) In addition to wing vertical bending and torsion 
flexibility, fore and aft bending flexibility was included. 

(2) The inboard nacelles were considered as dis- 
tributed masses on flexible struts. Displacement due 
to strut extension was considered negligible. 

(3) The outboard engine was considered as a dis- 
tributed mass rigidly attached to the wing. 

(4) The stabilizer was considered rigid. 

(5) Flexibility of the body was included. 


Aerodynamic 

(1) Gusts were considered as two dimensional and 
symmetrical. 

(2) Gradual penetration of the airplane into the gust 
front was considered in computing wing and stabilizer 
air forces. 

(3) Kussner and Wagner lift growth functions were 
used in computing air forces. 

(4) The effects of finite span on the circulatory part 
of the air forces was accounted for by a modified Weis- 
singer analysis. 


TABLE 1 
Computed Symmetrical Airplane Vibration Modes 


Mode Frequency, 
Number cps Description of Primary Motions 

1 1.3 First wing bending 

2 2.2 Inboard nacelle strut side bending 

3 2.8 First wing fore and aft bending and 
nacelle strut vertical bending—some 

4 3.0 wing torsion 

5 3.6 First body vertical bending 

6 5.6 Second wing vertical bending 

7 6.0 Nacelle strut torsion 

8 oF First wing torsion and second body ver- 


tical bending—some second wing ver- 
tical bending and first wing fore and 
aft bending 

9 7.2 First wing torsion and second wing ver- 
tical bending—some second body ver- 
tical bending and first wing fore and 
aft bending 
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Fic. 2. Typical experimental gust record 

(5) The stabilizer was idealized as a single lifting 
surface element, and the effect of wing downwash was 
accounted for by modifying the stabilizer angle of 
attack, 


a, = [1 — (de/da) la, 


(6) In computing the generalized aerodynamic forces, 
it was assumed that both the Kussner and Wagner lift 
growth functions for all stations on the wing and sta- 
bilizer could be expressed in terms of distance traveled 
with respect to a single reference chord instead of the 
local chord, as follows: 


WV (Vt,/b)) = VA Vt,/b,), KA Vt/b) = K(Vt/b,) 


The reference chord was selected to minimize the total 
lift growth error for the Kussner Function. 

(7) It was originally assumed that there were no air 
forces acting on the fuselage. In light of early results, 
which showed the theoretical short period pitching fre- 
quency to be too high, it was necessary to revise this 
assumption implicitly by modifying the aerodynamic 
coefficients that correspond to the rigid body motions. 
The coefficients were corrected to agree with rigid 
model wind tunnel data. 


In early studies, undertaken to investigate the effects 
of speed, altitude, and gross weight, an analog com- 
puter was used. Because the large size of the problem 
taxed the computer capacity, it was necessary to com- 
pute loads by the approximate, mode displacement 
technique (described in detail in reference 6) which 
required substantially less equipment than the direct 
approach of summing air load and inertia load com- 
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ponents (the load summation method). While the re- 
sults obtained with this technique were adequate for 
the trend comparisons, its accuracy was considered 
questionable with regard to providing data to compare 
with experimental results. Technique comparison 
studies were made on the digital computer which sub 
stantiated this belief. All subsequent analysis utilized 
digital computers and the load summation process. 

The theoretical analyses for the specific gust com- 
parisons were made using elastic modes 1, 4, 5, and 6. 
These analyses were for an airplane weight of 125,000 
Ibs., for a c.g. location of 26 per cent mean aerody- 
namic chord and for sea level altitude. Experimental 
values varied about +2 per cent of the weight and about 
+5 per cent of the mean aerodynamic chord on c.g. 
location from those used in the analyses. These vari- 
ations were considered to have only minor effects on 
the results. The average pressure altitude during 
the gust course runs was about 2,400 ft. This discrep- 
ancy was also considered to be of secondary impor- 
tance. The airplane forward speed and the gust 
length were obtained from the gust record. Loads were 
computed for the stations shown in Fig. 3. These were 
not the same ones for which the experimental loads 
had been obtained. In the comparison of theoretical 
and experimental results, therefore, the theoretical 
data were plotted and interpolated for values at the 
experimental stations. 


RESULTS AND DISCUSSION 


Technique Evaluation 

Initial theoretical studies were undertaken to check 
the effects of various assumptions and techniques. 
Analyses were made using, in addition to rigid transla- 
tion and pitch, mode 1 only, modes 1, 2, 3, and 4, and 
modes 1, 4, 5, and 6. The shears, bending moments, 
and torsions were computed for 1 — cos shaped gusts of 
various lengths by both mode displacement and load 
summation methods. At least one similar investiga- 
tion has been made of the B-47 airplane. Results are 
reported in reference 6. However, the assumptions 
made there appeared to be sufficiently diiferent to 
warrant the undertaking of the present study. 

Figs. 4 and 5 show present analysis results typical 
of inboard and outboard stations. The peak incre- 
mental bending moments plotted are as functions of 
gust length for a forward speed of 300 knots. The 
weight c.g. and altitude are the same as for the specific 
gust comparison analyses. Results for modes 1, 2, 3, 
and 4 are not shown because they were virtually the 
same as for mode 1 only. At the inboard station the 
maximum peak moment predicted by mode displace- 
ment analysis using mode 1 only is about 12 per cent 
less than that predicted using the load summation 
method. A similar comparison for modes 1, 4, 5, and 
6 shows mode displacement values about 8 per cent 
below the load summation moment. The peak values 
predicted by the load summation method using mode 1 
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only appear to be only about 10 per cent below those 
of the larger analyses. 

At the outboard wing station the variation due to 
technique is more substantial. It can be seen that, for 
both the mode 1 only and the mode 1, 4, 5, and 6 analy- 
ses, the Joad summation process results in maximum 
peak loads which are about 31 per cent greater than 
those predicted by the mode displacement method. It 
is interesting to note that the results of the analysis 
including the higher modes shows a pronounced peak 
at about 10 chords, indicating a critical phasing of the 
gust with one of these higher modes. The effect of 
this higher mode is apparently overestimated by the 
mode displacement method, as evidenced by the fact 
that the moment predicted drops off very rapidly, 
after its early peak, to a value substantially below the 
mode J] only mode displacement value at a gust length 
of 50 chords. For gust lengths greater than about 20 
chords the load summation analyses results for the mode 
1 only and the modes 1, 4, 5, and 6 are again in good 
agreement, the latter giving from about 8 to 15 per 
cent greater moments. 

As mentioned previously, loads computed from 
modes 1, 2, 3, and 4 analyses were essentially the same 
An examination of the modes provides 
Mode 2 is the airplane 


as mode 1 only. 
an immediate explanation. 
mode which principally involves side bending motion 
of the inboard nacelle strut, and mode 3 is that mode 
which is predominantly wing fore and aft bending. In 
both of these modes, the generalized force due to the 
gust is relatively small. The response buildup requires 
several cycles, therefore. Consequently, the modal 
peaks occur sometime after the peak loads are passed. 


Specific Discrete Gust Comparisons 
A careful review of the low level gust data resulted 
in the selection of six gust encounters as being suitable 
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Fic. 6. Comparison of experimental maneuver and gust incre- 
mental wing shear. 
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Fic. 7. Comparison of experimental maneuver and gust incre 
mental wing bending moment 


for theoretical comparison. The fact that the gust 
shapes were not as idealized in the analyses was con- 
sidered relatively unimportant on the basis of results 
presented in reference 3. 

In order to determine the magnitude of the dynamic 
effect, experimental maneuver load per unit c.g. accel- 
eration values were determined from data obtained 
previously on the same airplane for symmetrical roller 
coaster maneuvers. So that the effects shown would be 
principally dynamic, these data were obtained for the 
identical gross weights, speeds, altitudes, and airplane 
c.g. locations as the experimental gust data. To facili- 
tate comparison, the experimental gust incremental 
loads per unit c.g. acceleration were divided by corre- 
sponding maneuver values. This ratio is, of course, 
commonly referred to as the magnification factor. 
Moreover, it is evident that the results are a direct 
indication of the potential accuracy of the steady-state 
theory in conjunction with the familiar gust formula, 


AN = K,U,Vm/(W/S) 


where  & = an alleviation factor 
U, = gust velocity 
V = airplane speed 
m = airplane lift curve slope 
W/S = wing loading 


The results are presented in Figs. 6 and 7. The 
comparison of the shears shows remarkable agreement 
between the gust and maneuver values at the two 
stations for which shears were measured. It can be 
seen that the root average magnification factor is less 
than one, and that the trend is toward an increasing 
factor with span. On the other hand, the bending 
moment magnification factor plot shows a substan- 
tially larger dynamic effect. It can be concluded that 
the shear magnification factor at outboard stations 
would also be much greater than one. The effect of 
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wing bending moment. 


speed and gust gradient variations are clearly illus- 
trated by the large data scatter. The use of a steady- 
state theory to predict gust loads would obviously re- 
sult in substantially underestimating the magnitude 
of the bending moments over a large portion of the 
span. 

The trend toward increasing magnification factor 
with span should not be surprising. It is pointed out 
in reference 7 that this result is obtained for transient 
loading on a simple unrestrained beam. It is worth- 


while to note, however, that theoretical B-47 analyses 
show a leveling off of magnification factors outboard of 
the outboard nacelle (W.S. 694.4). 
in fact, there is a decrease at the outboard stations. 
Comparison results of theoretical analyses are pre- 


In some cases, 
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sented in Figs. 8 and 9 in the same form as the maneuver 


load comparison. As mentioned previously, loads 
were computed with the load summation method using 
elastic modes 1, 4, 5, and 6, as well as rigid body pitch 
and translation. With regard to the shear ratios, the 
theoretical values at station 120 are, on the average, 
less than the experimental values, and the trend is a de- 
This, it 
will be remembered, is opposite to the trend shown by 


creasing ratio with increasing wing station. 


maneuver data comparisons. 

An inspection of the bending moment ratios shows a 
substantial improvement over those based on steady- 
state theory. It can be seen that the dynamic analysis 
results are in reasonable agreement with the experi- 
mental data over the entire portion of span surveyed. 
Inboard of station 400, the average ratio lies only 
slightly below one. At station 505, however, this aver- 
age decreases to about 0.95, and it rises to about 1.05 
at station 642. 

A possible explanation for some of the analysis dis- 
crepancies can be seen by examining typical load time 
histories. Figs. 10 and 11 are plots of experimental and 
theoretical incremental bending moments at wing sta- 
tions 120 and 642. It can be seen that both experi- 
mental traces show evidence of higher mode contribu- 
tions not shown by the theoretical results, the effect 
at the outboard station being more pronounced. The 
effect of the higher mode can also be seen in the shear 
loads, even near the root, as shown by Fig. 12. From 
the data shown, it would appear that this higher mode 
contribution might be as large as +5 per cent, depend- 
ing upon the spanwise station. This would account 
for the discrepancies noticed in the average values. 

The frequency of the higher mode is about 7.2 cps. 
Referring to Table 1, it can be seen that there are two 
wing modes in this frequency range. It is possible to 
surmise, then, that, in order to obtain better agreement 
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THEORETICAL AND 
with the experimental data, one would have to include 
modes 8 and 9 in the analyses, bringing the total to two 
rigid body modes and six elastic modes. Since the 
seventh mode involves very little motion other than 
torsion of the nacelle strut, its exclusion would result 
in essentially no error in the wing loads. 

The reader will note that two sets of points are miss- 
ing from the shear comparison plots, Figs. 6 and 8. 
The original data were mislaid prior to reduction. 
However on the basis of the corresponding bending 
moment data, one might expect no marked difference 
in comparative results. 


Gust Velocity Results 


Both gust probe and gust vane data were reduced to 
provide a measure of the gust velocities for the specific 
gusts analyzed. Theoretical gust velocities were com- 
puted from the c.g. acceleration increment for the 
appropriate conditions by the simple ratio 


U, theory — ANr AN, 


where AN; is the peak c.g. acceleration from the flight 
test records for a specific gust and AN, is the peak c.g. 
acceleration from the analysis for unit gust velocity and 
the experimental conditions. The results are shown in 
Fig. 13 as a plot of theoretical gust velocity versus 
experimental gust velocities computed from both the 
gust vane and the gust probe data. It can be seen 
that both the probe and the vane values are generally 
higher than those predicted by the theory. The un- 
corrected probe data showed a large velocity component 
due to boom bending. When this component was re- 
moved by curve smoothing, the results were brought 
into better general agreement with the vane data. 
The discrepancies experienced are probably due in 
part to neglecting pitching effects in computing the 
gust velocities from the experimental angle of attack 
data. However, it is possible that these discrepancies 
are also due in part to the fact that the theory may 
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mental bending moment time histories at wing station 642. 
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predict c.g. accelerations for a given gust which are 
greater than those actually encountered. Considering 
the fact that the peak loads per unit acceleration 
agreement is good, a prediction of larger c.g. acceler- 
ations than actually occur would result in a conserva- 
tive prediction of the gust loads. 
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CONCLUSIONS 


(1) The steady-state theory is inadequate to accu- 
rately predict discrete gust loads on the B-47. Loads 
predicted may be less than 50 per cent of the actual 
loads for some spanwise stations. 

(2) Using presently available methods, the peak 
B-47 wing loads due to gust encounters can be predicted 
accurately by dynamic analyses. 

(3) The general shapes of the theoretical load time 
history curves are in good agreement with the experi- 
mental results. However, the experimental results 
show a substantial higher mode contribution not evi- 
dent in the theory. 

(4) The theory, as applied in this study, may pre- 
dict c.g. accelerations somewhat larger than actually 
occur. 

(5) The loads computed by the mode displacement 
method are below the experimental loads for the range 
of conditions covered in the present study. This dis- 
crepancy is largest at the outboard stations. 

(6) The nacelle strut side bending mode and wing 
fore and aft bending mode due to discrete gust en- 
counters have a negligible effect on the wing vertical 
loads. 


APPENDIX——-A METHOD OF THEORETICAL ANALYSIS OF 
THE Dynamic Gust LOADS PROBLEM 


General Considerations 


The equation of transient motion of the aircraft due 
to symmetrical vertical gust excitation is written in 
terms of the motion and external air load perturbations 
relative to the dynamic state of the airplane prior to 
encountering the gust. 

The continuous distribution of mass throughout the 
airplane is idealized by a finite number of concentrated 
masses located along the wing and fuselage elastic 
axes, as indicated by Fig. A-1. In the idealized dy- 
namic system, the 7th concentrated mass is assigned in- 
ertia properties that are identical to those for the cor- 
responding segment of the actual system when referred 
to the 7th reference point on the elastic axis. Struc- 
tural stiffness and elastic deformation of the airplane 
are represented by normal vibration mode frequencies 
and displacements at each inertia load reference point. 
At any particular instant of time, the net perturbation 
of all concentrated masses in the idealized system is 
expressed in terms of normal mode spatial data and the 
time-dependent generalized coordinates. 

The air forces acting on the airplane are expressed 
in terms of lifts along the quarter chord and pitching 
moments about axes perpendicular to the stream direc- 
tion at the quarter chord air load stations shown in 
Fig. A-2. 


Equation of Motion 


The equation of aircraft motion is derivable from 
fundamental work and energy theory. For symmetri- 
cal vertical gust excitation, m elastic and two rigid 
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Fic. A-1. Airplane elastic axis system and inertia reference 
stations. 


body degrees of freedom are arranged in the following 
order to describe the transient motion of the airplane: 


General- 


ized 
Degree of Modal Coordi- 
Freedom Description Value nate 
(1) First symmetrical 
normal vibration 
mode pi? g: 
(n) nth symmetrical 
normal vibration 
mode oi Yn 
(n+ 1) Rigid body vertical 
translation 1p @tyt Gn4t 
(xn + 2) Rigid body pitching 
rotation | pnt} Gn+2 


In terms of the above modal data and generalized co- 
ordinates, the displacement perturbations of the con- 
centrated masses in the idealized system are given by 


tel Fe 
iss = [6] 195 (1) 
where all linear displacements are nondimensionalized 
by dividing by a reference length 5,. All angular dis- 
placements are in radians. 
In equation of airplane motion in terms of dimen- 


sionless time is 
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1(V/b,)? [6’Jo] {q"}| + 1(V/b,)? X 
[o’ Joi] [(wis/wr)*] tg} = (OF (2) 
where dimensionless time is defined by 


S= Vt/b,, wo, =V/), 


[o’J@| and [@’J¢i;| [(wii/w,)?] are the generalized in- 
ertia and stiffness matrices, respectively. 

The left-hand side of Eq. (2), when set equal to zero, 
is the equation of free undamped vibratory motion of 
the airplane in the elastic degrees of freedom, in addi- 
tion to unaccelerated vertical translation and pitching 
rotation of the center of mass of the system. The right- 
hand side of Eq. (2) consists of the generalized air forces 
which are jointly dependent on the gust and subsequent 
aircraft motion. 

The generalized air forces }Q} are derived from the 
expression for virtual work done on the airplane by 
the air forces acting at the wing and horizontal stabilizer 


air load stations. 


Unsteady Air Forces Due to the Gust 

From incompressible thin airfoil theory, the lift on a 
streamwise strip of an infinite wing due to a gust of 
any arbitrary shape may be expressed in the following 


manner : 
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and a, = U,(S)/V (5) 
For tapered plan forms, the reference semichord /, 
is selected so that the error in lift growth over the entire 
wing is a minimum. It then follows that the Kussner 
lift growth function approximated in reference 8 reduces 
to 


W(S) = 1 — 0.5e~°-'85 — 0.5e-4 (6) 


For tlie (1 — cos) gust 


{1 — cos 2r(b c) (S — s,)] (7) 


U,(S) = (1/2) U,... 
where s; is the time lag between the instant of gust front 
penetration by the first and 7th air load panels, as shown 
by Fig. A-3. 

— &)/b, (S) 


(S —s,) =O for 02> (S — s;) > c/b, (9) 


V) sin 2r(b,/c) (S — s;) (10) 


Unsteady Air Forces Due to Aircraft Motion 

The unsteady lift and moment due to arbitrary 
motion of a two-dimensional wing based on incom- 
pressible thin airfoil theory are given in references 9 
and 10. For a streamwise strip of width AY, the un- 
steady lift and moment are 


bl = 2rpV*b,3(AY/b,) [(b/b,)K*h” + 
(b/b,)?K*a"” + (b/b,) K*a’| + 
mp V*b,*(AY/b,) [(b/b,)*h” + 
(1/2) (b/b,)'a” + (b/b,)?a’} (11) 


m = mrpVb,2(AY/b,) [(1/2) (b/b,)*h” + 
(b/b,)'a’ + (3/8) (b/b,)4a”| (12) 


where 
Ki = [ K(S — a)h"(a)do (13) 
K*a’ = { K(S — a)a"(a)do (14) 
K*a’' = f K(S — a)a’(a)da (15) 


The Wagner lift growth function, approximated in 
reference 8, is also written in terms of the same refer- 
ence semichord used in Eq. (6), 


K(S) = 1 — 0.165e~%-°%5 — 0.335e—-5 (16) 


General Unsteady Air Forces 


The total unsteady lift and moment on the stream- 
wise strip due to the gust and transient motion are ob- 


tained by combining Eqs. (3), (4), (11), (12), (13), 
(14), and (15). 
Soll = pV%b,5 ([Pi] fv*a,’| + [Ps] JK*h"| + 
(m | 10 f \K*a"f 
[P;] j0 (+ [Ps] Jh"l + [Ps] JOU (17) 
\K*a’f la” la’ : 
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where 
AY 2 0 b 
7 0 ¥ a (18) 
b, b, 
AY 
0 0 0 0 O 
b, 
Py] = A} 0 2r Ol] b (7) (19) 
[2] = b, b, b, 
9 AFii0 ojjo 0 
b, 
a A} 27 O b 
b, b 
AY 
0 0 0 oe @¢ 
b, 
(7) l (Py . 
b, 2 NO; =) 


AY 0 b\? 
[Ps] = 0 ||" 0 (;) (22) 
b, b, 
O b\? 
F (7) 


Prior to using this equation in the theoretical analysis 
it will be modified to include the effects of airfoil thick- 
ness and camber, compressibility, sweep, and spanwise 
variztion of section lift coefficient. 

The theoretical value of the lift curve slope, 27, in the 
circulatory lift terms is replaced by the two-dimensional 
zero Mach Number lift curve slope ad. The effects of 
compressibility and sweep for the given airfoil and plan 
form will be included in the circulatory lift terms by 
multiplying by an equivalent Glauert correction G, 
which is defined below. Following the work of Gray 


b, 
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Fic. A-3. Initial position of the airplane relative to the gust 


front. 


and Schenk," the effects of finite span, plan-form taper, 
and sweep on the spanwise lift distribution are inserted 
in the circulatory lift terms by introducing a finite 
span induction matrix [.S;~'] into Eqs. (18), (19), and 
(20). Derivation of the elements of the induction 
matrix is based on a modified Weissinger lifting line 
theory and is described in detail in reference 11. 

Rewriting Eq. (17) in a form to include the above 
modifications, and all lifting surface elements on one 
side of the wing plan form, 


fol| = C, -_ he S¥* a9" + Gow [Pale SK*h"\ + GeolPale JOU + [Ps] fh" + [Ps] JO L) (23) 
lm Fs 10 f K*a"f \K*a’l la" la’f 

where Co = pV*b,5 (24) 
Gem = 1/V1 — M? 08? Ayayy (25) 

AY r T° 
ae ={[2°]fo]]][s--][o] || [J [o] | [2] 25 
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[S,] = [6/2] [Si] (31) 


where [.S,] is the downwash matrix described in refer- 
ence 11. 

Eqs. (23) and (25) through (31) are also applicable to 
the stabilizer when the subscript w is replaced with s. 
The effects of unsteady downwash on the stabilizer are 
omitted and the angle of attack at the horizontal tail 
is expressed in terms of that for the wing by 


a, = [1 — (de/da) la, (32) 
C.= pVb,?[1 — (de/da) | (33) 


The vertical displacement and pitching rotation of 
the air load stations along the wing and stabilizer 
quarter chord may be written in terms of the normal 
mode spatial data and generalized coordinates. 


Sh\ 


+ ee FZ ]w i) (34) 
ji a 
c = [H].{q} (35) 


The use of the subscripts w and s to indicate wing 
and stabilizer air loads is clear. Therefore the sub- 
scripts will be discontinued in equations which are 
applicable to either the wing or stabilizer. 

After substituting Eqs. (34) or (35) into Eq. (23), the 
lift and moment in terms of the generalized coordinates 


are 
fol = C(G[Ui] Je* al + G.[U2] | K*q"} + 
im § 10 f 


Ge[Us] {K*q’} + [Us] fg"! + [Us] 'a'1) (36) 


[Ui] = [Pi] 
[U2] = [Po] [77] (37) 
[Us] = [Ps] | 


Generalized Forces 
The work done on the airplane by the air forces 
during a virtual displacement is 


6bW = —2[6q] [Hu Jbl — 2[6q] [H'], Jott (38) 
im 4. lm 4, 
The generalized forces indicated by the right-hand side 
of Eq. (2) are 
LQ} = {Qu} + {Qs} = ((0/d8q)6W} (39) 
After combining Eqs. (34), (36), (37), (38), and (39) 
the generalized air forces are 
{Q} = —2C (G, [Ri] f¥* a, + 
10 f 
G.[Re] {K*q"} + Ge[Rs] | K*q’} + [Rs] ig") + 


[Rs] 'a'1) (40) 


where 
[Ri] = [H’] [Pi] 


[R2] = [H"] [P2] [7] (41) 


[Rs] = (H'] [Ps] a) 


Equation of Aircraft Motion Due to a (1-cos) Gust 
The equation of aircraft motion in dimensionless 
form is obtained by combining Eqs. (2), (39), and (40) 
and dividing through by the scalar 2p V’*),’. 
[Ne] {9°} + [Mil ta’} + [No] is 
[Re] (K*q") + | 
[Ri le jy*a,’ ( 


10 ‘B Yo 
where 
[No] = (1,/2b,>) [6'J¢] + [Ralw + 
[1 — (de/da)| [R4|, 
[Mi] = [Rsle + [1 — (de/da)] [Rs] ; 
[No] - I, 2pb,°) lo’ Joi: | [(wi: w )?] 5 (43) 
[Ro] = Gey [Relw + Ges[1 — (de/da)] [Re], 


| 
[R;] _ Gew[R3 lu + G, {1 — (de da) | [R 
[R:] = Gew [Ri Ju 
[Ri], = G..{1 — (de/da)] [Ri] 


Transient Structural Loads 


The transient shears, bending moments, and torsions 
due to the gust and disturbed motion of the airplane 
may be computed when the elastic deformation of the 
structure isknown. The ‘‘mode displacement” method 
described in reference 7, relates the internal forces at a 
given location in the aircraft structure with the general- 
ized displacements in the elastic degrees of freedom. 
Shears, bending moments, and torsions in the wing 
at the locations shown by Fig. 3 are obtained from the 
following equation: 
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g{S)i i= 1,2,...,2 (44) 


M(S) 


jig 
tris) §, \. 


While the mode displacement method gives only a 
rough approximation of the transient structural loads 
when used with a small number of elastic degrees of 
freedom, a minimum of computing equipment and time 
is required. 

A more precise method of computing the transient 
shears, bending moments, and torsions in a continuous 
system, while using only a finite number of elastic 
number of elastic degrees of freedom, consists of adding 
the effects of inertia loads and external air loads, which 
Using the “load summa- 
bending 


(vis)| = \ 
| 


Ss 
b 
t 


are computed separately. 
tion’ method, the net incremental shears, 
moments, and torsions on the wing at the locations 
shown by Fig. 3 are obtained from the following 
equation: 
VS) | = 
M(S) 
lms) }., 


Computation of transient structural loads by the 
method requires approximately the 


” 


(S)} + [D]ufbl\ (45) 


\m f w 


[C lu 1g” 


” 


“load summation 
same amount of computing time as the solution of the 
equation of motion when the problem is solved by digital 
computation. When operational amplifier analog com- 
puting techniques are used, solution of the load equa- 
tion by the “load summation” method requires ap- 
proximately the same amount of equipment as required 


” y j , Ao uv 
1g pre i — [Vi] Va + es Qn \ 


n : 1 c 
[Val D K(Sy — 0) fa'}do — 5 [Val dan! 
i=0 ~ 


where 
[Vi] [Q-"] [Mi] 
[V2] [Q-*] [No] 
[Vs] [(Q-"] [Re] . 
[Ve] = [0-*] [Ri] sd 
[Velo = [(0-*] [ile 
[Vs], = [Q] [Ai]. , 
[Q] = [No] + (Ao/2) [Mi] + (Ao /2)? [No] + 
(Ao/2) [Ro] + (Ao/2)? [Rs] (51) 


For the case where the airplane is flying straight and 
level prior to gust entry, and the gust vertical velocity 
profile is of (1 — cos) shape, the initial conditions for the 
generalized coordinates are 


q"| =0 for S<0 


f -. we 

tg = ta} (52) 

The solution of the equation of motion [Eq. (42) ] for 
the time histories of the motion perturbations is ob- 


tained by solving the finite difference equation, Eq. 


SPACE 


A 
— [V2] a. + Acq,’ + (5 
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for solution of the equation of motion. A highly ac- 
curate computer is required for this method of load 
computation since the inertia and air load contributions 
to the net incremental loads are numerically large and, 
in general, are of opposite sign. 

The acceleration of any point on the airplane may be 
computed from 


AN(S) = [®] }q¢"(S)} (46) 


A Finite Difference Solution of the Equation of Motion 


The desired solution of the equation of aircraft mo- 
tion [Eq. (42)] consists of the time histories of all per- 
turbations of the generalized coordinates over a time 
interval of sufficient length to include the peak dynamic 
response of the airplane to gust excitation. 

A finite difference method for solving the equation of 
airplane motion by digital computation is outlined 
below. 

Numerical integration of the functions g”(S) and 
q'(.S) over the interval nAc < S < (n + 1)Ao, using the 
trapezoidal rule, gives the following finite difference 
equations: 


(47) 
(48) 


Q'ntt = Qn’ + (Ao/2) (qn” + Gg" n41) 
= gn + Acq,’ 2)? (a," rel 


where (” + 1) refers to the (7 + 1)th assigned value of 
S for which the equation of motion has just been solved 
for g”. After substituting Eqs. (47) and (48) into 
Eq. (42) and collecting terms, the equation of motion 


Qn+1 


becomes 


o 2 p 
) an’ aot Ag — 


7 K(S,, _ a;) 1q”} 


i=0 


o . \ - n 7 j a am ) 
n A _ V; w (3; we i) Pa in 
q ng Ae [Vs] ay 7) Yo f a 
n fa’) 
5|s (S, — o;) : Ao (49 
Wy ee 1 
(49), for }q"(S)} at S = Ao, 2Ac,..., nda, (n + 1)Aa, 


. over the desired interval of time. After solving for 
{q”(S)} at each assigned value of S, {q’(S)} and |g(S)} 
7) and (48). 


Computation of the incremental structural loads and 


are computed from Eqs. (4 


airplane accelerations at each assigned value of S may 
be carried along as the solution of Eqs. (49), (47), and 
(48) continues. Where the mode displacement method 
is used, the structural loads and accelerations at each 
value of S are obtained immediately from Eqs. (44) and 
(46). 

Where the load summation method of computing the 
structural loads is used, it is necessary to compute the 
external air loads for each assigned value of S through- 
For S = (n + 1) the net lift 


and moment on each surface element are given by 


out the time history. 
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[U's] 50" | nt T [Us] fa'}ess) (53) 


Using the load summation method, the incremental 
structural loads on the wing at each assigned value of 
S may now be computed by substituting the results 
obtained from Eqs. (49), (47), (48), and (53) into 
Eq. (45). The airplane acceleration time histories are 
computed by the use of Eq. (46). 
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An Approximate Analysis of Nonlinear Flutter Problems 


(Continued from page 32) 


tered. More important, perhaps, is the advantage 
that, by reducing to an equivalent linear problem, the 
role of the noniinearity is clearly exhibited. A feeling 
of knowing why the system behaves the way it does is 
certainly prerequisite to an efficient attempt of making 


any improvement. 


REFERENCES 


1 Woolston, Donald S., Runyan, Harry L., and Byrdsong, 
Thomas A., Some Effects of System Nonlinearities in the Problem 
of Aircraft Flutter, NACA TN 3589, 1955. 

2 Woolston, Donald S., Runyan, Harry L., and 
Robert E., An Investigation of Effects of Certain Types of Struc- 
tural Nonlinearities on Wing and Control Surface Flutter, Journal 
of the Aeronautical Sciences, Vol. 24, No. 1, pp. 57-63, January, 
1957. 

3 Kryloff, N., and Bogoliuboff, N., Introduction to Nonlinear 
Vechanics (a free translation by S. Lefschetz), Princeton Uni- 
versity Press, Princeton, 1947. 

‘Fung, Y. C., The Flutter of Buckled Plate in a Supersonic 
Flow, GALCIT OSR TN 55-237, July, 1955. 

5 Hedgepeth, John M., Flutter of Rectangular Simply Supported 


Andrews, 


Panels at High Supersonic Speeds, 1AS Preprint No. 713, January, 
1957. 

6 Stoker, J. V., Nonlinear Vibrations in Mechanical and Elec- 
trical Systems, p. 115; Interscience Publishers, New York, 1950. 

7 Popov, E. P., On the Use of the Harmonic Linearization Method 
in the Automatic Control Theory, NACA TM 1406, January 1957. 

8’ Kochenburger, R. J., A Frequency Response Method of Analysz- 
ing and Synthesizing Contactor Servomechanisms, Trans. AIEE, 
Vol. 69, Part I, pp. 270-284, 1950. 

Johnson, E. C., Sinusoidal Analysis of Feedback-Control 
Systems Containing Nonlinear Elements, Trans. AIEE, Vol. 71, 
Part II, pp. 169-181, 1952. 

1 Theodorsen, T., and Garrick, I. E., Flutter Calculations in 
Three Degrees of Freedom, NACA Report No. 741, 1942 

1! Fung, Y. C., Flutter of Curved Plates With Edge Compression 
in Supersonic Flow, GALCIT OSR TN 57-187, March, 1957. 

12 Shen, S. F., and Hsu, C. C., Analytical Results of Certain 
Nonlinear Flutter Problems, Journal of the Aeronautical Sciences, 
Readers’ Forum, Vol. 25, No. 2, pp. 136, 187, February, 1958. 

13 Woolston, D. S., and Andrews, R. E., Remarks on “ Analyti- 
cal Results of Certain Nonlinear Flutter Problems,’’ Journal of the 
Aeronautical Sciences, Readers’ Forum, see p. 51, this issue. 

14 Shen, S. F., Author’s Reply (to reference 13), Journal of the 
Aeronautical Sciences, Readers’ Forum, see p. 52, this issue 





An Improvement of the Myklestad Method 
for Flexural-Vibration Problems 


S. MAHALINGAM* 
University of Ceylon 


SUMMARY 


The paper presents a method of obtaining the next approxi- 
mation for the natural frequency by means of a simple auxiliary 
calculation based on the results of a Myklestad! table. The 
basis of the method is the change of natural frequency caused by 
a small change in either the mass or the elasticity of one com- 
ponent of a flexural system with a finite number of degrees of 
This improvement results in a rapid convergence of 
conditions are con- 


freedom. 


successive approximations. Various end 


sidered. An illustrative example is given. 
SYMBOLS 
a = amplitude of slope 
I = second moment of area about flexural axis 
1 = length of beam segment 
m = mass of particle 
M = amplitude of bending moment 
Q = amplitude of shear force 
s = total number of beam segments 
T = maximum kinetic energy in the system 
V = maximum strain energy in the system 
x = distance measured along beam 
y = amplitude of deflection 
w = a natural circular frequency of the system. The sub- 


scripts 1 and 2 attached tow denote the first and second 
approximation, respectively. 


INTRODUCTION 


N DEALING with problems of free vibrations of non- 
I uniform beams the Rayleigh-Ritz and Stodola’s 
methods provide satisfactory solutions for the funda- 
mental frequency. However, these methods become 
very cumbersome when the higher natural frequencies 
have to be determined. For instance, in order to ob- 
tain the mth natural frequency, the first (7 — 1) nat- 
ural frequencies have also to be calculated. A direct 
numerical method using a Holzer type of iteration was 
first given by Myklestad in 1944, and soon afterwards 
a similar method, differing only in detail, was given by 
Prohl.? Since then an alternative method has been 
given by Bellin’ while a simplification of the Myklestad 
method was given by Bishop‘ recently. But the essen- 
tial characteristics of all these methods and modifica- 
tions are the same. The nonuniform beam is replaced 
by a system consisting of a set of point masses con- 
nected by sections of light elastic beams. A frequency 
is assumed, and, starting from one end of the system, 
the beam equations are satisfied at each station. If 
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the assumed frequency is approximate, then one of 
the boundary conditions at the other end will not be 
satisfied. The procedure then reduces to constructing 
the table for several trial frequencies and determining 
that frequency for which the boundary condition is 
In general, this involves a great deal of 
numerical computation. 

An improvement of the Holzer method which pro- 
duces a rapid convergence of successive approximations 
was given by the author’ recently. In this method, 
the next approximation for the natural frequency is 
given by a simple calculation based on the results of the 


satisfied. 


previous Holzer table. 

A similar improvement of the Myklestad method is 
given in the present paper. Flexural vibration prob- 
lems are inherently more complicated than torsional 
ones since changes of slope and deflection are involved 
and additional difficulties are caused by different 
boundary conditions. Consequently, the next approxi- 
mation for the natural frequency is obtained by con- 
sidering a small (fictitious) change either in the mass 
of one particle or in the elasticity of one segment of the 
beam, depending on the boundary conditions. 


EFFECT OF A SMALL CHANGE IN A FLEXURALLY 
VIBRATING SYSTEM 


Change of Mass 


Consider a system of point masses m, ms, . . . Ms+1 
connected by light prismatic beam segments of lengths 
li, lb, . .. 1, and second moments of area J, Jo, .. . J, as 
shown in Fig. | 

When the system is vibrating freely at a natural fre- 
quency w, let V and 7 represent the maximum strain 
and kinetic energies, respectively. Then V = 7. 

The maximum strain energy is that stored in the 
beam when it is acted upon by static loads P;, Po, ... 


P.+, where P; = mw*y;. If m, is a variable mass, then 
V = SM, Yo, - - » Vets) 
T = F(y, Yo, . . - Vet1, w?, M,) (1) 
= w*Fi(y1, Vo, . . - Vsti, M,) j 





ser ”, ™,; M7, be 
O00 ——__» x 
(s+) s 3 2 2 2 f 
iso & ad p> 4+ 4 4 
J 
Fic. 1. Flexural system with a finite number of degrees of 


freedom. 
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We then have 
s+1 s+1 
- > (OV /dy,éy, = >> Pi by, = 
t=1 t=1 


s+1 
> mw*ydy; (2) 
1 l 


s+1 
Since r= (1/2) > mM wy; 
#=1 
OT/dy, = Mwy; ) 
OT /O(m,) = (1/2)w*y,? (3) 
OT /O0(w?) = O(w*F;)/O(w?) = T ot 
and 
s+ 
67 = >> (OT /dy,)by; + [O07 /0(m,) ]6(m,) + 
i=1 
[OT /O(w?) | 6(w?) 


ed} 
= DY) mwy oy; + (1/2)w29,2 6(m,) + 
i=1 
(T/w*)6(w?) 
Equating 6V and 67 we get 


5(w*) = —o"} [6(m,)y,? 1/¥ m yt} (5) 


This result is similar to that obtained when there is 
a small change in the moment of inertia of one rotor 
of a torsional system® 


Change of Elasticity 


If J, is variable, we have 


V = f'(m, Ve, Vst1, I) | 
T = F'(y, y2 - Vst1, w*) 7 (6) 
= w* F,'("1, V2, Vs I 
6V = >) mwy by; + (Of'/O/,)6/, (7) 
i=1 
5 +1 
6T = > mw’y by; + (T/w?)d(w?) (8) 
t=1 


Equating Eqs. (7) and (8) we get 


s+1 
av[[a /2) p> miy | (9) 


AV = (Of’/ol,)él, 


6(w?) = 


where 


AV represents the change in the strain energy because 
of a change in J,, the deflections being kept constant. 
An alternative expression for AV may be obtained as 
follows: 

Consider a set of point loads W;, We, . Ws+1 pro- 
Then the strain 


ducing deflections y;, yo, . . . Vs+1. 


energy is 


st1 
2) > Wa: (10) 

i=1 
This can also be written as V = $(Wi, Wo, .. . We+1 
i) 

When /, is given a small change, let 6W;, 6W2, 
5IV,+, be the changes necessary in the loads to keep the 
deflections the same. The change in strain energy is 


then 


s+1 . 
AV = DY (0¢6/0W)sW, + (d¢ o1,)81,) 
i=1 


s+1 (11) 
= ¥ yéW, + (0¢/0/,)4/, ( 
*=1 
But from Eq. (10) we have 
+1 
AV = (1/2) > 6W, y, (12) 
i=1 
It follows from Eqs. (10) and (12) that 
AV = —(0¢/0/,)é/, (13) 


Since ¢ is a function of the loads and J/,, 0¢/0/, can be 
obtained easily. In the expression for ¢, J, occurs only 
in the term representing the strain energy of the rth 
segment. The variation of bending moment between 
stations in linear so that 


V, = strain energy in rth segment of the beam 


lr 
= (1 2E1,) f M*dx 
0 
= (] 2E1) [ [M, + (x/l,) (Mii — M,)]? dx 
0 
= (Il,/6EI,) [M,? + M,M,+1 + M,+;7] 
from which 0¢/O0l, = —(V,/I,) 


and substituting in Eq. (9) 


s+ 
6(w?) = vat, |a 2) _ my? (14) 
/ i=1 


is the strain energy in the 
If, however, we consider a 


In the above expression V, 
rth segment of the beam. 
more complex system like that of a nonuniform beam on 
flexible supports where the stiffness of one of the sup- 
ports is variable, then V, would be the strain energy 
in that support. 

Eqs. (5) and (14) form the basis of the improvements 


outlined below. 
IMPROVEMENT OF THE MYKLESTAD METHOD 


Free-Free Beam 


Consider a free-free system A consisting of masses 
s + 1) as shown in Fig. 1, vibrating at 


mle = 1,2... 
Then, at any station ” we have 


a natural frequency w. 


n 
0, = dm; wy; 
i=1 


n—1 


> m; wy; (x; 


#=1 


I 


M,, — Xn) 

If the slope a and deflection y; at the station | are 
assumed, then the slope and deflection at any station 
n may be expressed as linear functions of a and 4; as 


i, = Q1Dn — Vidn’ t 
Fa * — ayn + Win’ J 


The functions ¢,, ¢,’, ¥,, and y,,’ are connected by the 


relations 


(15) 
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Pnt+1 _ Pn =o k,Bn ie &,,'3,,’ j 
nti. = Gn’ + Rie i SG, (16) 
, , ) 
Wrt1 = Vn + LnOnt+1 i 2nB» = Be B,, ( 
Wnt1’ = v." + bi onti’ = 2nGn aig Ca Ge 
where 
Ri = ba GE ie,’ - be Fer. £n = 
LV/Ghle, £4 = 17/241, 
and 


n Nn 
B, = Do mwy:, Gr = > mw,’ 
t=1 i=1 


n—1 n—1 


B,! = Lomw ix; — xn.) = DB: 
Ppa i=1 


n— 1! of 
Gn’ = > m, wy;'(x; — x2) = LG, 

i=1 i=1 
The notation used above is that given by Timoshenko.® 
When the Myklestad table is constructed for the system, 
the boundary conditions require zero bending moment 
and shear force at the station (s + 1). 

For the first condition we have 


M +1 = —aB,+,' + WGst1' = 0 


which gives 


an/V1 = Get’ / Bot’ (17) 
The deflections may then be written as 
Yn = Ml—(Get1'/Bst in + vn] 
giving 
Q, = p> mw" [bi’ — (Go+1'/Bs+1') Wi] 
= wilGn — (Got1'/Bst+1')By] 
and Osta = ¥1[Get1 — (Get1'/Bs+1')By+1] (18) 


If the Myklestad table is constructed using an ap- 
proximate frequency w, then Q,+; will not be zero. 
However, this boundary condition can be satisfied by 
replacing m,+, by a suitable mass m,+,’.. Then w will 
be the exact natural frequency of the modified system 


B for which m,’ = m,fors = 1, 2,...s and 
Ms+1' = - (Q, O17? Vs+1) 
The change of mass is therefore 
Smt = Met — Msi! = Qyt1/anr2yct1 (19) 


The system B can now be converted to the original 
system A by adding 6m,+, to the last mass. There will 


Vs4 1” 
Wo41 


Va" + le bots” — BDs — Bs'D,’ 
Vs + 1s bst1 — g0B, — ge'B, 
Vs" +e bs” + (lsksDe + lsks'Ds’ — Ds — 8.'D.') 


1939 
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then be a corresponding change in frequency 6(w,") giv- 
ing the second approximation 
wo? = wy? + 6(w)”) 


ae 
= w? — ws'(sm. Vet v/ 2 mi'ye) 
1 I 


; 
/ 
' | (> P ; 
= wi? — | QstiVsts mii? — 6ms+Ws+1? 
i=1 (20) 
st+l i 
= |(= mM 4"y i? — 2Q0.+19s4 ) = 
s=1 
s+] 
(= mM, on?V i? — OstiVs4 | 
i=1 


The calculation of w.* from Eq. (20) is quite straight- 
forward and similar to that given by the author for} 


7 eater preg 


torsional vibrations. 


Simply Supported Beam 


If R; and a are the support reaction and slope, re- | 
spectively, at the station 1, then the slopes and de-| 
flections at other stations may be expressed as linear | 
functions of these two quantities as 


an= AiPn = Rid,” t 
7s = an — Rw," 


The functions ¢,, ¢,", ¥», and y,” are connected by 
the relations 


(21) } 


Ont1 = On + RnBn + kn'Br’ 
uti” = bn” + haDn + by'Dy’ wid 
Wnt = Vn + Ln@n+1 = £rB,, — ga'B,’ ge 
Wnt 1” — Yn” + lndn+1” a £rDp, a g.'D,’ 


where B, and B,’ have the same meanings as before, 


and 


D, =1+ > mwy;" 
+=] 


n—1l 


D,’ = LD, 
1=1 


At the station (s + 1) the beam is simply supported and 
Vst1 = sti — Rids+i” = 0 

a) we (23) | 

The bending moment at this point is 


My+1 = Ri[Deti’ — (sti Ys, 


Qa! R, _ Y; 


giving 


)Byt1'] (24) 





If the assumed frequency is approximate then 
M+, ¥ 0, but this condition can be satisfied by making 
a small change in J,. 


We have from Eq. (22) 


¥. + 1.6, + (2B, + 1,k,'B.’ — Bes — g. B,') 


In order to satisfy I/,+,; = 
I, reg | =e BI,’. 
Then w, is the exact natural frequency of the modified system B, and from Eq. (24) 


0, let 7, be replaced by /,’ 


so that J,/J,’ = (1 + 8) where 8 is small, giving, 
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MYKLESTAD METHOD 








PIV. 


20) 





2) 


re, 





giving B 


By+,' 


FOR FLEXURAL-VIBRATION PROBLEMS 


Devi! _ Wa" + oda” + (1 +B) (lskwDs + Iyks'D,! — gD, ~ g,'D,) 





¥. thos + (1 + B) (kB, + 1k,'B,’ — 2,B, — g,’B,’) 





= sts” + B(,k.Ds + I,ks'D,' —_ gD, — g;'D,') 
Vsti + BRB, + 1ks'Bs’ — g.B, — gs'B;') 


The modified system can be converted into the orig- 
inal system A by adding @/,’ to the second moment 
The corre- 
sponding change of natural frequency given by Eq. 
(14) is 


of area of the sth segment of the beam. 


s+] 


5(w;?) = wav,’ /( 2) > mwn’y;”? 
i=1 


amplitude of displacement and V,’ = strain 
The second ap- 


where y’ = 
energy of sth segment, of system B. 
proximation is therefore 


s+1 
we” = wi} (© mM wy ;'? + 23V,) 
i=1 


Beam on Flexible Supports 


s+l1 
9 19 
> M 4iW1"°Vi ~ 


i=1 


(26) 


For an accurate determination of the critical speeds 
of steam turbine rotors the flexibility of the supports 
must be taken into account. The forced vibrations 
of such a system have been considered by Miller.’ 

Let c; and c,+,; be the stiffnesses of the supports. If 
the slope and deflection at the station 1 are assumed, 
then those at other stations can be defined by relations 
similar to Eq. (15). Also 


n 
> Mwy; 


i=1 


Qn 


— ay 


n—1 n—1 
YS mwy; (xi — Xn) — Mn DA 
i=1 


+=! 


M,, 


I 


The functions ¢,, ¢,’, Wn, and y,’ will be connected 
by the relationships in Eqs. (16), with B, and 5 
having the same meanings as before and 


n 
G, = > mw,’ — « 


i=1 


n—1 


Gc." = » LG; 
i=1 


The boundary conditions are M,+; = 0 and Q,+: = 
Cs+1Vs+1. When the assumed frequency is approximate, 
these conditions can be satisfied either by a change in 
m,+, or by a change in c,+;. If the latter is adopted 


Vi! = (1/2) 654+1'¥s+1?. 


Cantilever Beam 


Using the free end as station 1, the slopes and de- 
flections at all stations are expressed in terms of a; and 


' y, and the boundary conditions at the other end are 





satisfied by a change in /;. 
If the fixed end is chosen as station 1, then the slopes 
and deflections at all points are expressed in terms of 





Vsti M +1 


~ Ry[Byt1'(leksDs + lgke’'Ds’ — geDz — B0'Ds') — Det1'(leksB, + leks'B,’ — g,B, — g,'B,')] 


R; and M,, and the boundary conditions at the other 
end are satisfied by a change in m,+. 


Beam With a Sliding End 


The concept of a beam with a sliding end is very use- 
ful in the study of frames. A member in a structure 
may be so constrained by its adjacent members that 
one end or a point in it may move as if it had a slider 
attachment. 
vibrating symmetrically, the midpoint of the hori- 
An 
analysis of beams having this end condition was made 
by Bishop* 
When such a beam has a uniform cross section, the re- 
ceptances are calculated by direct methods while non- 
uniform beams require the use of the Myklestad method. 


For instance, when a portal frame is 


zontal member may be regarded as a sliding end. 


9 


in a study of the vibration of frames. 


In the present discussion, if the station (s + 1) is at 
a sliding end, the boundary conditions can be satisfied 
by a change in m,,). 


GENERAL REMARKS 


it is seen from the above examples that, in order to 
satisfy the boundary conditions, a change is made in 
that component of the system which is the last to 
enter the beam equations. Problems involving a 
moving end (free, sliding, or flexibly supported) re- 
quire a change in the last mass. A fixed or a pinned 
end requires a change in the elasticity of the last seg- 
ment of the beam. The auxiliary calculation is simpler 
in the case of a change of mass, and this is to be pre- 
ferred when the choice is available. 


NUMERICAL EXAMPLE 
The numerical example given by Timoshenko’ will be 
considered here. To determine the natural 
frequency of a uniform beam of negligible weight carry- 
ing three equal masses, as shown in Fig. 2, 


second 


m = 1.5 Ib. in.—! sec.? 


EI 
l; =,=~);,=-),= 


= 40 X 106 Ib. in.? 
10 in. 


For this mode of vibration station 3 is a node, and the 
exact natural frequency can be calculated by consider- 
ing a simply supported beam of span 20 in. carrying a 
central mass m. This method gives w» = 400. In 
order to illustrate the application of the present method, 
an approximate value of w;? = 150,000 will be assumed. 


The Myklestad table will then be as follows: 
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Myklestad Table 


n mw? XK 10-6 Wn B, x 10°* 5,’ X 10" Pn+1 

l 0 0 0 0 1.0000 
2 0.225 10.0000 2.2500 0 3.8125 
3 0.225 29.3750 8.8594 22.5000 20.5117 
4 0.225 132.5391 38.6807 111.0988 117.2016* 
5 38.6807 497 . 9004 


0 718.7503* 


Using Eq. (25) we get 8 = 0.27016. This gives E/,’ 
= [(40 X 10®)/1.27016] = 39.49213 X 10° Ib. in.? 

In Table 1, the quantities marked with an asterisk 
are those affected by the change in J, and the figures 
given refer to the modified system. Since M,+,' = 0, 


V,! = (44/6EI,') (M,’?) 
= 5.02605 X 10-6 Ib. in. 
s+1 /s+1 
ww? = wf (X myw,"y,'? + 2av.') 2 meosty | 
i=1 i=1 


150,000 X [(34.66912 + 2.71566) K 10~§ + 
34.66912 * 10-6] 


161,750 


It is thus seen that from a first approximation having 
an error of about 6 per cent in w’, and using only one 
operation of the Myklestad table, we get a second 
approximation which has an error of about 1 per cent. 
In the normal procedure an error of this order would 
have been obtained after about two or three tables. 
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37 .2396 
178.6621 
920 .5494* 
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f 
id 
vn” X 108 D, D,,’ dni” X 10% y,'/Ri X Ws M,'/Rif 
1.0000 0 1. 2500 0 ' 
4.1667 1.93875 10.0000 6.1719 8.64097 f 
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50.5154 132.5391 148.7007* —8.91089 —9.74518 
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5 a 3 Pg / 
f~ - = —( )— 
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and Hsu! relates to the flutter of 

of freedom 


— ARTICLE by Shen 
two and three degree 
nonlinear spring characteristics 
analytical procedure whereby the structural nonlinearities may be 
and treated by 


systems with certain 


Their article concerns an 
equivalent linear stiffnesses 
The procedure is applied to 


approximated by 
usual methods of flutter analysis 
some of the examples which we have treated by analog computer 
methods.?? 4 

Shen and Hsu raise serious questions concerning the 
of our results in that they find different flutter boundaries from 
those which we have indicated. An appraisal of the situation 
indicates that two different approaches are and that 
The so-called discrepancies 


accuracy 


involved 
our published results are correct 
noted in reference 1 are due simply to the fact that reference 1 
compares two dissimilar quantities. 

The two quantities involved are illustrated in Fig. 1 which 
represents a typical analog computer record of the limited ampli 
tude flutter of a nonlinear bending-torsion system. The quan 
tity ap (which we employ as a flutter indicator) represents an 
initial displacement from equilibrium, applied at t = 0. The 
quantity apg (employed by Shen and Hsu) represents the limit 
amplitude of the resultant response. The quantities a9 and a; 
are not directly related and cannot be compared. 

The approach of Shen and Hsu assumes the existence of an 
established flutter oscillation of prescribed amplitude a, in the 
degree of freedom which contains the nonlinearity and seeks the 
at which the specified oscillation will be 
are concerned with the response 


critical velocity sus- 
tained. We, on the other hand, 
of the system to a prescribed, abrupt displacement, ao, representa 
a sudden control movement. At a 


tive of a gust encounter or 
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Fic. 2. Variation of preload, 5p, and flutter amplitude in tor- 


sion, ar, with velocity ratio. 
reference 2. 


Bending-torsion example of Fig. 4, 


given velocity, we seek the magnitude of the disturbance required 
to induce flutter. Our analog approach will give, as well, the 
magnitude of a; of a resultant flutter of limited amplitude. 
At a given velocity, the values of a» obtained by the two ap- 
proaches should agree. 

In reference 1 certain incorrect assumptions are made concern- 
This seems due, in part, 
to an incomplete reading of our text and, in part, to the fact 
that we did not define the systems fully. Fig. 1 of reference 1, 
for example, treats our bending-torsion system with free play and 
preload (location of the equilibrium point relative to the flat spot) 
and regards preload as constant. In reference 2 we tried to in- 
dicate the preload was variable by the statement, “‘. . . the preload 
was comparable in effect to a deflected tab and produced a 
moment which varied with the velocity.”” In order that proper 
comparisons may be made, we present in Fig. 2 the variation in 
preload with velocity employed in our computations (ordinate 
at left of figure). 

As a result of the questions raised in reference 1, we have 
redone’ our computations for the bending-torsion 
example and been able to reproduce our published 
results. In addition, we have made some observations of the 
limit amplitude, a,, in the more violent of the two flutter modes. 
In order that flutter comparisons will be possible, the variation 
with velocity of the amplitude of a, is included in Fig. 2 (ordinate 


ing the properties of our examples. 


analog 
have 
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at right). 
“experimental”’ data point presented by Shen and Hsu at ap = 
2.5° in their Fig. 1. Actually, this is a symbol from the legend 
of our figure (see Fig. 4 of reference 2) and should be deleted 

In reference 1, computations are also made for a case corre- 
sponding to our example of wing-control surface flutter with free 
When 
proper attention is given to the distinction between an initial 


A more easily resolved point of issue concerns an 


play and preload in the aileron degree of freedom.’ 


displacement 8) and a resultant amplitude of oscillation By, 
and, when the same amount of preload is used (we used dp = | 
0.7 67), the results of the two approaches should be in substantial f 
In reference 3, we included the amplitude of the 
flutter oscillation of the control surface on Fig. 6. It is this result 
Incidentally, 


agreement. 


with which comparison should have been made. 
reexamination of our results shows that we should have presented 
an amplitude ratio of Bp/éy = 0.7 at V/Viin = 1.0 rather than a 
value of zero. 

Both the analytical approach and the analog computer ap- 
proach to the handling of nonlinear systems appear equally valid } 
provides significant information. The method of |} 
Shen and Hsu predicts the existence of a flutter oscillation and its | 


and each 
amplitude. The analog approach provides this information and, 
in addition, predicts the magnitude of disturbances which the 
system can sustain before a flutter oscillation develops. 
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F* PHYSICAL PROBLEMS, careful experiment is always the 
ultimate check mathematical 
accurate numerical evaluation serves to confirm or disprove 


of theory. For solutions, 
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We cannot agree more with Woolston 
“Remarks” that 
computer and our analytical procedure! are complementary 


approximate analysis. 


and Andrews in the above their analog 


to each other. Indeed, we rather hoped that our predictions 
would fall in line with the analog results. 
ing to us that these authors apparently support our position 
that our method might be successul in calculating flutter bound- 


If is therefore hearten- 


aries of nonlinear systems. 

We of course apologize for having inadvertently taken their 
parameter in the flutter boundary plots as the steady oscillation 
amplitude. With the newly published Fig. 2 of the note above, 
where 6p increases from about 0.01° for V/Viin & 0.35 to about 
0.08 for V/Viin & 0.8, the comparison with our theory in the 
form of ag vs. V/Viin is now shown as Fig. 1, replacing the 
The agree- 
In the other example, for which 


previous comparison given as Fig. 1 of reference 1. 
ment seems quite reasonable. 
they referred to the lower right of Fig. 6 of reference 2, we are 
unable to establish the coordinate for 8,/6 after correcting the 
zero point (6f/6 = 0.7 at V/Viin = 1.0), as suggested in the 


above note. Thus no comparison is possible. 
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The Effect of Slip on Induced Pressures} 


J. A. Laurmann* 
Research Physicist, University of California, Berkeley, Calif. 
July 29,1958 


Ss“ TSIEN’S HISTORIC PAPER of 1946! on superaerodynamics, 
writers concerned with low Reynolds Number effects have 
Tsien’s classification of free-molecule, 
transition, slip, and continuum flows. But for those writing in 
the field of rarefied gas dynamics, one of the most frustrating re- 
sults has been the continuing failure to find an experimental 


used the concepts in 


situation of the external-flow type in which slip effects are, in 
fact, observed.2 This stituation exists partly because of the re 
duced mean free path length behind the shock associated with 
supersonic flow past a body, and partly because of the masking 
of possible slip effects by the interaction of the surface boundary- 
layer flow with the bow shock wave 

Recently, however, measurements of the induced pressures on 
a flat plate* at high values of the hypersonic interaction param 
eter and low Reynolds Numbers have shown for the first time an 
effect that could be attributed to slip. 
experimental points obtained at two Mach Numbers in the low- 
density wind tunnel at the University of California for the ratio 
of the induced pressure p to the free-stream pressure p., versus 


Thus, Fig. 1 shows the 


the hypersonic interaction parameters x = M*/+/Re:, where M 
is the free-stream Mach Number and Re, the local plate Reynolds 
Number based on free-stream values. Every interaction theory 
so far developed for high values of x predicts an approximately 
linear dependence of pressure on x. The fall of the experimental 
points below a straight line is what would be expected from slip 
at the wall, as is also the earlier drop-off for the smaller Mach 
Number. 


+ Abstracted from The Momentum Integral Method Applied to the Com- 
pressible Boundary Layer with a Slip Boundary Condition, University of 
California, Project Rep. HE-150-160, 1958. Work conducted under spon- 
sorship of the Office of Naval Research and the Office of Scientific Research 
at the University of California, Berkeley, Calif. 

* Now at the Space Technology Laboratories 
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In order to check this possibility, a slip analysis was made of 
the zero pressure gradient boundary-layer flow, using a Pohl 
hausen type solution, combined with the Stewartson Transforma 
tion to give the compressibility effect. The approximate solution 
thus obtained predicts qualitatively, though not quantitatively, 
the behavior described above. A principal difficulty has been to 
take correct account of the interaction between the shock wave 
and the boundary layer, and, for checking the validity of the slip 
assumption, this difficulty was circumvented by using a small 
perturbation analysis to obtain the change in induced pressure 
caused by slip. Since, as mentioned before, all the no-slip high 
x interaction theories predict an essentially linear dependence 
of pressure on x, the experimental results have been compared 
with slip theory by subtracting the slip correction to the pressure 
from an empirical no-slip straight line obtained from the experi- 
mental points. Fig. 2 shows the results plotted, not with the 
customary variable x as abscissa, but versus M?/y ‘Rez, a param- 
M)/[(p — po)/pPao] exclusively de- 
pends in the hypersonic case, M— ©. The result of the revised 
method of plotting is to better correlate the results for different 
Mach Numbers, and it can be seen that the slip theory gives at 
least the qualitatively correct variation of induced pressure with 
Mach and Reynolds Numbers. The point of departure of the 
slip-flow theoretical curve from the no-slip straight line agrees 
namely, M?/+/ Rez > 


eter on which the ordinate (1 


with the criterion given by Sherman‘ 
0.4 in this case—for the breakdown of continuum flow 
check on the slip assumption cannot be made, however, until a 
more rigorous interaction theory has been developed, or other 


A closer 


experiments, say in the subsonic range, have been performed. 
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On the Stable Shape of an Ablating Graphite 
Body* 


Daphne Christensen and Rolf D. Buhler 

Research Scientist and Project nee Director, Respectively, 
Plasmadyne Corp., Santa Ana, Calif 

August 7, 1958 


A“ INTERESTING PHENOMENON on the ablation of graphite in 
hot air has been observed. 

In preliminary testing in the plasma-jet hyperthermal tunnel, 
earbon rods have been observed to ablate to sharply pointed, 
cone-like shapes rather than to blunt-nosed bodies, as might have 
been anticipated. The pointed shape appears to be a stable 
geometry and no significant changes in contour were noted as 
ablation continued, as shown in Fig. 1. 
experiments were made in a 
having the following approximate properties: 
Mach Number 2.4, stagnation enthalpy 5,000 B.t.u./lb., free- 
stream temperature** 6,500°R., velocity 9,400 ft./sec., 
0.00021 Ib. /ft.*, Reynolds 


These one-in.-diameter, arc- 


heated air jet 


static 
pressure 0.5 psia, density Number 
3,300 per in. 

The 1/4-in. cylindrical carbon rod was placed coaxially in the 
free jet of the axisymmetric nozzle. Surface temperatures along 
the carbon rod, measured with an optical pyrometer, were of the 
order of 4,500°R but were not sufficiently accurate to establish 
reliable heat-transfer information. However, the weight-length 
correlations and sharp point formation were consistent for vari- 
ous test procedures, as shown in Fig. 2 

Similar results have been observed using argon and nitrogen 
as working fluids, although the ablation rate is so low that a com- 


plete ablation sequence, as shown in Fig. 1, has not been obtained. 


* This work was supported by the Army Ballistic Missile Agency, Hunts 
ville, Ala 
** Based on assumed thermodynamic equilibrium, 


for Air. 
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(a) PHOTOGRAPH OF ROD TIP BEFORE AND 
AFTER ABLATION 








(b)CARBON ROD ABLATION SEQUENCE AT 
ONE MINUTE INTERVALS 


Fic. 1. Ablation of a carbon rod in air. 





‘SPACE 


SCIENCES—JANUARY, 1959 
























































12 06 os 7 
ee | 
°o BUA 5 j 
10 10c00 BTU/, 3 —~ 
e | — — 
77] 
G) 2 | 
Be] 3. = t t 4 
z A ; 4 
z WEIGHT 7 a 
5 a ~ ss 
te. cn ™ | 
4 2 — a 
i | oe 4 
oe 
2 Ol —-— af 
i 
a 
ft) oO 
0 1 2 3 4 5 6 7 


TIME IN MINUTES 
Fic. 2. Ablation of two graphite rods in air. 
Nevertheless, sharply pointed conical shapes produced with air 
remained sharply pointed while ablating in either the argon or 
Tentative measurements indicate that the weight 
ft., was about 0.10 with air 


nitrogen jet. 
loss rate of graphite, in lb. /min.-sq. 
as compared to 0.02 with nitrogen and 0.004 with argon. 
Experimental and theoretical research on sublimating bodies 
is continuing. The measured ablation rates have been compared 
with theoretically calculated values based on the laminar flow re- 
lations with chemical reactions, as given by several investigators 
In these calculations, the wall shear stress has been taken such 
as to fit the approximately known carbon wall temperatures, 
whereas the cold-wall (Blasius) shear values have generally been 
other authors. Satisfactory agreement between these 
ablation rates and the 
Attempts are also being made to explain theoretically 


used by 
calculated 
obtained. 


the unusual equilibrium shape described abeve 
/ 


measured values has been 
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Note on the Shock-Induced Unsteady Laminar | 


Boundary Layer on a Semi-Infinite Flat Plate* 


S. H. Lam and L. Crocco 

Research Associate, Gas Dynamics Laboratory, and Robert H. 
Goddard Professor, Respectively, Princeton University, 
Princeton, N.J 

July 28, 1958 

L’ T A SEMI-INFINITE FLAT PLATE be immersed in a quiescent 

gas of infinite extent and let its leading edge be the origin 

of a coordinate system (x, y) with the plate occupying the posi- 

tive x-axis. A plane shock wave perpendicular to the x-axis, 

advancing in the positive x direction with constant velocity U,, 

For 


reaches the leading edge of the plate at the instant ¢ = 0. 
t > 0, the shock travels over the plate. A boundary layer will 
develop on the plate between the leading edge and the moving 
shock. The gas velocity behind the shock is denoted by um, and 
the strength of the shock is characterized by the number A which 
is the ratio of the shock velocity U’, to uo 

The purpose of this note is to study this unsteady boundary 
layer under the following assumptions: (1) the boundary layer 
(2) the boundary-layer approximations are 
(4) pw = constant, and (5) 


is laminar, valid, 
(3) the plane shock remains plane, 
P, = uc,/k = constant. 

It is interesting to note that in the limiting case of a very weak 
shock, when A = o, the problem is mathematically equivalent 
to the impulsive motion from rest of a semi-infinite flat plate 


* This research is part of a program sponsored by the Air Force Office of 


Scientific Research, under Air Force Contract Number AF 18(600)-498. 
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Also, when the flat plate in question is 


treated by Stewartson.! 
infinite in both directions, ‘‘similar solutions’’ exist. This class 


of ‘‘shock-tube” boundary layers has been treated by Mirels,? 
Bromberg,*® Bershader and Allport,‘ and others 
In the present study, the Crocco transformation will be 


adopted. One defines: 


T = potty? [u(Ou/Oy)| (la) 

H = (h — hh)/ho (1b) 
a = x/Uolt (le 

B= u/u (1d) 

Y = Uo*t/v (le) 


7 


where subscript 0 denotes conditions in the free stream. The 
shear function + is assumed to be non-negative in the region of 
interest. The compressible laminar boundary-layer equations 
for zero pressure gradient can be written as: 
Momentum equation 
r’rg8 = ty + (1/7) (8 — a)ta (2: 

Energy equation 
72[uo?/ho + (1/Pr)Hge] = 

H, + (1/7) (8 — a)Ha + (1/2) (7?)gHell — (1/Pr)] (2b) 


The region of interest is a rectangular right prism in the 
(a, B, y) space open in the + y direction: 
Oo< et 4A, OSs tl, vy 20 (3 


The pertinent boundary conditions for the shear function r+ 


and the enthalpy function H are: 


initial condition y = 0, +r— ©, H bounded (4a) 
leading edge a=0, r—~o, H bounded (4b 
shock a=A, r—~o, H bounded (4c) 
wall 8 =0, 7g = 0, H or rHg given (4d 
free stream B=1, r=0, rrp =0, H=O (4e) 


It is noted here that boundary conditions exist at two a@ sta- 
the velocity and enthalpy profiles are required to be uni- 
Eqs. (2a) 


tions: 
form at both the leading edge and at the shock. 
and (2b) are singular parabolic differential equations.> Defining 
o(a, B, y) = T\/ay one obtains: 

¢*op8 + (1/2)B = alB — alga + ayo, (5a 


7 [(mo?/ho) + (1/Pr)Hpp] = a(8 — a)Ha + 
ayH, + (1/2) (¢)gHp [1 — (1/Pr)] (5b) 


The shock-tube ‘‘similar solutions”’ are: 


org, = 


"m= V/ (A — 1) y(A — a) M(B; A), 


—/a(A — 1)/(A — a) M(6; A) (6a 


1 
H, = —a(a) f, [11(0; A)/M(p; A)]!—~?" dp + 


f 


1 p ’ 
Pr( Uo? mw) J, ap f [M(A; A) M(p; A)]!~?* da (6b 
f 0 


where a(a) = Hyg(a, 0, y), and M(8; A) is the solution to the 
ordinary differential equation: 

M Mpg + (1/2) (A — B)/(A — 1) = 0 (7) 
with boundary conditions Mg(0; A) = 0, M(1; A) = 0.  Rigor- 


ous uniqueness theorems for each of Eqs. (2a) and (2b) are given 
by Lam.® It suffices to state here that for any given A > 1.0, 
this pair of similar solutions 7; and H, are unique in the S-region 
defined by: 


»: i<xtes 4, OS 8 <1, vy oO (8) 


which means that the presence of the leading edge at a = 0 intro- 
duces no perturbation in S and that all the effects of the leading 
edge are felt only in the I-region defined by: 


I: 0O<a<l1, 0<B6<1, y>0 (9) 
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Fic. 1 Numerical solution of ¢(a@8) for A = ©, 


The function 1/(8; A) has been integrated numerically for 
several discrete values of A > 1.0 by Bershader and Allport‘ and 
Bromberg.* It is noted® that the special case of A = 0 corre- 
sponds to the Crocco solution of the steady flat plate, 

M(B; A = 0) = o(a = O, B, y) (10 


A useful property® of the functions 1/(8; A) is that for any A > 
1.0, 

M(B; A) > M(g; 0 (11) 

One now has a well-formulated and correctly posed boundary 

value problem for ¢@ and H in the I-region. Ignoring the energy 

equation for the moment, the pertinent boundary conditions 


for @ in the I-region are: 


o(0, B, y) = M(B; 0) (12a 

¢ (1, 8, y) = M(B; A (12b) 

og (a, 0, y) = 0 (12¢c) 

o(a,1, y) = 0, ddg(a, 1, y) = O (12d) 
¢ (a, B, y) bounded (12e 


The term ¢, may be dropped from Eq. (5a) thus reducing it to a 
partial differential equation in only two independent variables 
aand pg. The boundary conditions to be satisfied by @ as listed 
(12) are “elliptic’”’ in character as a direct result of the 


in Eqs 
(5a).°® Equation (5a) can be 


singular parabolic nature of Eq 
written in the form [imposing boundary conditions (12¢c and d)]: 


' if) et” 
g(a, B) = Ql¢] = \f. ap | Iola, p)/d(a, A)| X 
B 0 


/ 1/2 
[A — 2a(A — a) (ha/) (a, d)|dr¢ (13 


which defines the integrodifferential operator Q. A sequence of 


iterants is now defined: 


o’(a, B) = 


max. of {.(8;0); Wa(A — 1)/(A — a) M(B; A)} (14a 


o"*) = {m(a)o™ + Olo™]}/[1 + m(a)] 14b) 
It is clear that if the iterants converge, the converged limit will 
(14) and will match 


The function 


be a solution to Eq. (5a) satisfying Eqs 
smoothly with the shear solution ¢, at a = 1.0 
(14b) is a weighting factor which was introduced to 


m(a) in Eq 
Preliminary calculations for the case 


stabilize the iterations. 
have been carried out on a Datatron digital computor 
The iterants tend slowly toward a reasonable 


A= « 
using m = 8a. 
function after about 12 iterations and fluctuate with an amplitude 
of less than one per cent until about the 18th iteration. The 
14th iterant is shown in Fig. 1. When further iterations are 
performed, the iterants appeared to diverge, especially in the 
region where a is close to one. Further work is being carried 
out to improve the iteration procedure to produce a more satis 


factory solution. 
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Considering the 14th iterant as a sufficiently accurate approxi- 
mation to the exact solution, it is seen from Fig. 1 that across the 
line a = £, the solution is well behaved. Across the line a = 1, 
the solution joins smoothly with the shock-tube similar solutions. 
It can be shown easily that the line a = 1 is a line of essential 
singularities in the sense that all a-derivatives are continuous 
across it, but the solutions on either side of it are not analytic 
continuations of each other. 

With ¢ known, the differential equation for H can be written as 


; E Pele, 0) PP 
H(a, B, y) = —a(a, y) dp + 
¢(a, p) 


B 
us f" , 
re f dp 
ho B 0 
f f° (\ — a)Hala, d, y) + yHy(a, », 7) 
aPr | dp 
B 0 


git Pr (a, r) go — Pra, p) 
where again a(a, y) = Hg(a, 0, y). 
h, is defined as the wall enthalpy when the wall is insulated, i.e., 
a(a, y) = O, the film coefficient is obtained immediately as: 


bla, dr) oa 


d(a, p) 


dX — 


dy (15) 


Since the recovery enthalpy 


f Cote (OT 
Se ee Oy Ju 
PolloC py " Os (16) 
2Pr f [p(a, 0)/d(a, y)]!- Pr dy 
0 
Cy = [p(OUu/Oy)] w/ porto? ] 


It is seen that when Pr = 1.0, the familiar result of Reynold’s 
analogy is again valid. The theoretical value of h, may be ob- 
tained by solving Eq. (15) by an iteration procedure similar to 
that employed for solving ¢. 
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On the Plateau and Peak Pressure of Regions of 
Pure Laminar and Fully Turbulent 
Separation in Two-Dimensional 
Supersonic Flow 


William J. Guman 

Research and Advanced Development Division, AVCO 
Manufacturing Corp., Lawrence, Mass., and Assistant 
Professor, Rensselaer Polytechnic Institute, Troy, N.Y. 


August 11, 1958 


hex PURPOSE of this note is to present simple closed-form 
semiempirical expressions for the plateau pressure ratio in a 
separated fully laminar two-dimensional supersonic flow and for 
the peak pressure ratio in a separated fully turbulent two-dimen- 


sional supersonic flow. The expression for the laminar case is 


arrived at by introducing the observation that the pressure ratio 
for separation is generally about one-half the plateau pressure 
ratio into an appropriate expression for the pressure ratio at sep- 
Analytically the final result reads 


aration. 


AERO/SPACE 
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Pp/bo = 1 + ClyM2/(M? — 1)'4*Rz/4 


where p,/po is the ratio of plateau pressure to the pressure at the 
beginning of interaction, J/ the free-stream Mach Number be- 
fore the interaction, and R, the Reynolds Number based on 
conditions before the interaction with the characteristic length 
L being measured from the leading edge to the disturbance causing 
The constant C is determined empirically from 
It was found that good agreement between 


separation. 
available data.! 
theory and experiment! was obtained by using a value of C equal 
to 1.1 for the cases of (a) a compression corner, laminar sepa- 
ration upstream of the corner but downstream of the leading 
edge, (b) a step corner in the order of several boundary-layer 
thicknesses, laminar separation upstream of the corner but 
downstream of the leading edge, (c) pure laminar separation 
induced by an incident shock. For a compression corner with 
laminar leading-edge separation, it was found that if a value of 
2.1 was used for the constant, that then predicted results checked 
well with available experimental data.'! In cases with transition 
in the separated layer the plateau pressure prior to the pressure 
rise at transition could be evaluated with fair to good correlation 
with experiment from the derived expression presented. 

Under the assumptions of simple-wave theory, the essentially 
wedge-shaped region of separated laminar flow is confined by an 
angle defined by the relation 


dy/dx = C (M? — 1)'/4/R,'/4 


where dy/dx represents the slope of the boundary confining the 
region of laminar separation and responsible for the pressure 
rise to plateau pressure. All other symbols correspond to those 
defined earlier. 

A simple semiempirical closed-form expression for determining 
approximately the peak pressure in a region of turbulent sepa- 
ration in a two-dimensional supersonic flow can be determined 
by a technique similar to that proposed for laminar separation. 
Using an expression for the separation pressure ratio in a turbu- 
lent flow? and an empirical relation between the separation pres- 
sure ratio and peak pressure ratio, one can derive the expression 

Pp/bo = 1 + (yM?/2b)(1 — K)/{1 + (vy — 1)/2]M?}) 
where p,/po is the peak pressure to the initial pressure, 17 the 
free-stream Mach Number ahead of the region of separation, k 
a constant approximately equal to 0.55? and 6 a constant roughly 
equal to 0.65. With these two constants, one can immediately 
evaluate the desired peak pressure ratio and compare the results 
with those of Fig. 3 of reference 2 or with results presented in 
reference 1. Fair to good agreement was obtained between 
predicted and measured pressure ratios up to a Mach Number of 
about 4.0. The flow deflection angle corresponding to the 
pressure rise in the fully turbulent case can be computed from 
the oblique shock equation relating the flow deflection angle to 
the pressure ratio across the shock and the Mach Number ahead 
of the shock. 
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The Oscillating Inboard Flap at Supersonic 
Speeds 


Edmund F. Tyler 
Engineer, North American Aviation, Inc., Downey, Calif. 
August 4, 1958 


— CASE of a control surface extending in from the wing tip 
along the wing’s trailing edge is treated by Miles! under the 
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restrictions: 

(1) Plan form of the control surface is rectangular. 

2) Free-stream Mach Number exceeds 1.0. 

(3) AR VM? — 1> 1(Mis the aspect ratio of flap). 

(4) The control surface is sealed at its edges. 

(5) Vibratory motion is independent of span station. The 
case of an inboard control surface will here be treated under the 
same conditions plus the condition that 


(6) Distance from control-surface tip to wing tip > control- 


surface chord / M2 — 1. 
For brevity Miles’ notation! will be used. 
over the wing where only the flap vibrates is given by 


—8x1(x, y, 2b AR — ye! (1) 


The vertical motion 


a(x, y,t) = 
for an angle of rotation @ of the flap. 
The corresponding expressions for downwash, potential func- 


tion, and pressure are 


downwash 
w(x’, y’!) = —(U/b) tan 6[(0/dx’) + tk’|z | 
= BU(1 + ik’x’)1(x’, y’, 2M — y’)\ 


(2) 


potential functicn 


= (bBU/x) 1(x’, x’ + y’, x’ + 2m — y’) X 


, 


f af dng(é, n) [1 + ik’(x’ — &)] (3) 
0 y'—2R 


ee 


where 
, > / 
g(f,n) = 1(E — |n}) [e—* Mg cos (Kk V/ — 9") Ve —) (4 
pressure 
p = —(pl tan @/b) [(0/Ox’) + tk’ ]¢ 
= ol 28 tan 0/m) 1(x’, x’ + y’, x’ + 2M — y’) X ‘ 
’ ny! (o) 
[ dtf(x’ — &) | dng(é, of 
/J0 J y’—-2R 
where 
f(x) = 5(x) + 2ik’ — kx (6) 


6(x) = delta function ( 


LIFT AND MOMENT 


The integral, c;, of the pressure 2M (omitting the constant 
—pU°8 tan 6) with respect to y’ from the left to the right leading 


Mach line is (see Fig. 1) 


2R+2x’ x 
cy = [1(x’) /2AMRr] ay’ | 
« i 0 


Interchanging the order of integration and integrating by parts 


gives 


G = [1(x’ ze) f déf(x’ — &) 
0 


a2 R+x’ 
| dng(é, n) + x’ 


* x 
2R +x’ \ 
f n dn[g(t, n) — g(t, — 2R)] ¢ (9) 


dng(t, n) — 


, 


—_— Fz 


Because of the step function 1( — |n|) in g(é, 7m), the first two 
integrals in brackets vanish while the third gives, putting g in 
explicitly, 


x 


Gy = [1(x’)/m] dt f(x’ — t)e~* ME x 
0 
é 7 ; 
f [dn cos («V2 — n2)/W# — n?] 


~§ 
The last integral of Eq. (10) is 


(10) 
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+ *- 
W ] = rr 
Fic. 1. Plan form of an inboard control surface 
g 
. dn cos (x’ “/ #? — n?*) 
(1/7) , = Jy (x't) (11) 
J -% Vi — 9? 
Whence Eg. (10) is 
(12) 


x 
Go = ue) f dt f(x’ — Bem ME Jo(x't) 
0 


It may be verified that Eq. (12) is the two-dimensional pres- 
sure. Therefore, the two-dimensional aerodynamic derivatives 
A-g and Aagg may be used to compute the lift and moment on 


a wing due to control surface motion 


HINGE MOMENT 


The integral c, of the pressure/2AR (except for a constant) 
with respect to y’ from the left to the right side of the control 


surface is, see Fig. 1, 


2R x’ 
[1(x’) ae) ay’ [ déf(x’ — —) X 
0 0 


2 
| dng(é, ) & (13) 
J y’—-2R 


Oo = 


Proceeding as before one obtains 


ax’ j 72 R 
Ch = 1(x’) } déf(x’ — &) ie 1) | dng(é, n) 
Fe 


22 R } 
(1 Rr) negli, n)dnpe (14) 
0 


Using restriction no. (3), Eq. (14) becomes, in view of the step 


function in g(£, ), 


ex 
or = ue) | dtf(x’ — t)e-™ ME y 
0 


: cos (x’ Vv ¢? — »*) ( 
dn ne 

/ > 9 

0 Ve — 9? 


The first integral in brackets was noted above; the second is 


(15) 


elementary, giving 


, 


x 
Ch = ue) f dé f(x’ — de 
0 


‘ Sp ; . , > 
| Jo(x’é) — (ARK) 'sin (x’E)} (16) 


ix’M¢t 4 


The first term of Eq. (16) as noted above is the two-dimensional 
pressure while the second represents a correction for aspect 
ratio. Noting the definitions of the two-dimensional aerody- 
namic derivatives in reference 2 in terms of the integrals 


1 
; = 
Sp = (QUti of un e~2*Mu J Oxen)du 
0 


it is seen that the correction for aspect ratio on the coefficient 
Apg for the moment of the control surface about its hinge is 


App = Apg® — [AR(x/2)V/ M2 — 1] 


(17) 


1 Apg* (18) 
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0 . ~ . . 
™ is the uncorrected coefficient and if 


1 
n ann One. . 
Sn* = (2", of un e~2*xMu cin (Qxu)du 
0 


then Apg* is the result of substituting s,* for s, in Apg 


where Apg 
(19) 


(OQ) 
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On the Solution of the Laminar 
Boundary-Layer Equations* 


Melvin Epstein 

Research Associate, Polytechnic Institute of Brooklyn, 
Brooklyn, N.Y 

August 4, 1958 


r HAS BEEN SHOWN! that the Prandtl boundary-layer equations 
can be written in the form 


f'' +af’+ea(l—f = 
(U/U.) 7 |f'(0f’/oe) — f"(Of/de)] (1) 
(the notation is that of reference 1), where 


a (Lg/U..) (d/dx) (Ug), 86 (L/U..)g2(dU/dx) (s 


bo 


* This research constitutes a portion of a general three-dimensional 
boundary-layer investigation supported by the USAF under Contract No. 
AF49(638)-217, under the guidance of Dr. Antonio Ferri. The author 
wishes to thank Dr. M. H. Bloom for his stimulating discussions of the 


problem 
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Fic. 2. Local shearing stress distribution on a circular cylinder 


For the purposes of this development, it is convenient to deter- 


mine the scale factor g such that a = 1. There results 


: 
g = (2U,/LL of U dx (3) 
0 


Eq. (1) can be cast into a nonlinear integral equation 
for the shear function f” by integrating between the limits 0 and 


n. Thus 


°7 
f’ — f'*® - | i(8 + Df? 4+ 
0 


(U/U x) g?|f'(Of’/oe) — f"(Of/2E)]} dyn — Bn — ff’ (A) 


Eq. (4) may be solved by iteration by introducing a suit- 
able approximation for /” in the right-hand side. 
expression represents a new approximation for f/” which may in 
turn be inserted into the right-hand side to yield a still better 
approximation to f”. The corresponding approximations to the 
velocity profiles and stream functions may be obtained by inte- 
gration. The value of f,”(0) is determined by requiring that 
f,”"( ©) = O (here n denotes the mth iteration). The closer the 
initial approximation is to the exact solution, the better will be 
each succeeding approximation. Hence, it is desirable that the 
initial approximation at least satisfy the correct boundary condi- 


tions of the problem. In addition, since Eqs. (1) and (4) are 


singular at infinity, the initial approximation should also have | 


the correct asymptotic behavior at infinity. The accuracy of the 
nth iteration can be improved if the initial approximation con- 
tains an arbitrary function of x which is determined by requiring 
that f,’(~) = 1. 
all of the boundary conditions of the problem. 

To illustrate the method, consider, for a first approximation, 
Substituting into Eq. (4) 


In this way, the mth approximation satisfies 


fy” = ae ® where a is a function of x. 


and requiring that f;"( ©) = 0 yields 


fi"(0) = [(1 + 38)/2a] — (U g?/2U..a?) (da/dé) (5) 
while a is determined from f;'( ©) = 1, i.e., 
[((5 + 7B) /4a?] — (5U g?/4U., a*) (da/dt) = 1 (6) 


In the 


case when similar solutions exist, da/d—é = 0, and the problem re- 


which may be solved by standard numerical methods. 


duces to algebraic manipulation. In particular, one obtains for 


fi”(0) 


38)/W5 + 7B (7) 


which is compared with the exact solution of the Falkner-Skan 


fi"(0) = (1 + 


equation in Fig. 1 and Table I. This simple analytical formula 
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is seen to give excellent agreement with the exact result for 
B> -(0.1. 

\s an example of the application of Eqs. (5) and (6) to a more 
general pré sblem, calculations were performed for a _ circular 
cylinder. The calculated variation of local skin-friction coeffi- 
cient is compared in Fig. 2 to the exact value calculated by the y 
The error at the position of maximum 


Blasius series method.! 
The calculated separation point is 


is less than 8 per cent. 
107° as compared to the exact value of 110°. 

The author has been unable to prove convergence of the iter- 
Successive iterations are straightforward but 
A second 


t 


ation procedure. 
become quite tedious, even for the second iteration. 
iteration was carried out for the case 8 = 0. The significant 
improvement (see Fig. 1) strongly suggests convergence. 

A somewhat better initial approximation has been investigated 
which is intermediate in labor to the first approximation dis- 
cussed above and a second iteration. In this case the initial 
is made to satisfy the differential equation at 


approximation 
© jin addition to the boundary conditions—i.e., 


» = Oand n = 


fy” = [(a? + B) (8) 


2a]le~—*" + [(a? — B)/2]ne—™ 


from which is obtained 
(11 + 358)a* — 28(1 + 58)a? — B%(1 + B) 
16a°® 


lla’ — 68a? 
16a® 


where a is obtained from 


fi"(0) = 


— 56? da 


dé 


Ug? 


a®>+ 8, 2 | 
U (9) 


8a> dé 


— B7(7 + 38) 


(73 + 109B)a* — 28(13 + 258)a? 
32a® 
Ug? | 73a4 — 46Ba? — 2782 da 5(a? + B) dg 
= : = + — =1 (10) 
l 32a? dé 16a® dé 


x 


The results of applying Eqs. (9) and (10) to the cases wherein 
similar solutions exist are shown in Fig. 1 and Table I. Excellent 
agreement is found for 6 > — 0.1. 

The method outlined above has been extended to the three- 
dimensional compressible laminar boundary-layer problem with 
A more complete description of the 


equally good results.” 


method will be presented at a later date. 
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. 7 Sat ee ee a oe 
| 
-a 
| Exact | f' = ae Eq. (8) | 
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| (0) T f ) ] Error f"(0) | Error } {''(0) T Error | 
- r i Sn a a a a ca 
| | | | | 
| | } | | 
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Comments on ‘‘Thermal Buckling of Clamped 
Cylindrical Shells’ 


David J. Johns 
Lecturer, College of Aeronautics, Cranfield, England 


August 15, 1958 


I* THE RECENT PAPER by Zuk,' an expression was presented for 
the critical buckling temperature of a clamped cylindrical 
shell in terms of the material and geometrical properties of the 
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shell. Restraint at the edges of the shell was assumed to be 


provided by rigid frames experiencing no temperature rise 

The assumption was made in reference 1 that the circumfer- 
ential stress induced in the shell when it experienced a temper- 
ature rise, 7, may be approximated by the function, 


oy = —(EaT/2) [1 + cos (27X/L)] (1) 


In other words, there is a compressive circumferential stress 
along the entire length, L, of the shell 
It is well known, that 
introduced at the junction of a shell and a rigid frame (or bulk- 
head) are extremely localized, and the circumferential stresses 
induced in the shell decrease rapidly away from the joint. In 
fact it may be deduced from reference 2, for similar conditions, 


that, 


however, the discontinuity stresses 


—EaTe~™* [sin mX + cos mX] 


~ 


where m = [3(1 — pw?)/r*t?]"/4 (3) 


From Eq. (2) it can be seen that a, becomes zero when mX 
equals 32/4 and then has a small tensile value for mX < (77/4) 
For all practical purposes however, the distance X (= 32/4m) 
defines the range of applicability of Eq. (1) and in Eq. (1) the 
symbol L should be taken to signify 327/2m rather than the length 
of the shell, unless L < 32/2m 

With this modification, the subsequent analysis of reference 1 
may be adopted yielding Eq. (5) of reference 1 for the critical 
buckling temperature. 

It should also be noted that the results of reference 1 can be 
extended to include the effects of flexibility of the restraining 
frames (or bulkheads). 

From reference 2 it can be shown that, 

o, = —[EaTe ™*/(1 + Rr,g)] [sin mX + cos mX] (4) 


where Ry is the ratio of shell and frame (or bulkhead) flexibili- 


ties 
For frame restraint 
Rp = 2(E(t/m)/ErA yr] (5) 
and for bulkhead restraint 
Rp = 2 E(t/m)/Egr|d/(1 Bp)| (6) 


where suffixes F and B refer to frame and bulkhead, respectively, 
A, is frame cross-sectional area, and d is thickness of bulkhead 
To summarize: the analysis of reference 1 has been modified 
to allow a more realistic form of the circumferential stress in the 
It has been found possible to include the effects of frame 
Hence the critical 


shell. 

(or bulkhead) flexibility into the analysis. 

buckling temperature is given by, 

Ter = |1 + Re.p] [w'!r2t(6Ln8 + 32L,,° 
{12(1 - pe?) ?*r2aL mtr tiw* 


2) + 384L,,{v] + 
— 16L,,7A?)] (7) 


where Lm = 3x/2m (8) 


and Rr.z is given by Eqs. (5) or (6). Eq. (7) may be solved in the 


manner suggested in reference 1. 
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A Note on Conical Flow 


Ernest C. Kennedy 

Senior Research Engineer, Convair, A Division of General 
Dynamics Corp., Daingerfield, Tex. 

August 18, 1958 


O- OF THE STANDARD REFERENCES dealing with conical flow 
is the well-known MIT Technical Report No. 1, Tables of 
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Supersonic Flow Around Cones, where the differential equa- 
tion of conical flow is solved numerically using the specific- 
heat ratio y = 1.405. The results are tabulated for various 
values of the conical ray angle. In almost all other publications 
dealing with supersonic flow of air, the value y = 1.4 is used. 

At this facility we recently calculated the flow fields around a 
series of biconic bodies of revolution using the method of character- 
istics. In the conical regime we used the MIT values (the only ones 
available), but in the characteristic region we took y = 1.4. Thus 
in the solution of our problem we employed two different values 
of y, 1.405 and 1.4. In many other reports that I have read, 
the same thing has happened. 
generally make the remark, “‘This slight difference in y does not 
affect the final results appreciably.” 

The object of this note is to see just how much 
is. In the tables below, us, a measure of the velocity of flow 
along a cone element, is taken as the independent variable. The 
other parameters are: My = free-stream Mach Number, 0, = 
conical shock angle, Mp, and 6, = Mach Number and flow 
direction immediately behind the shock wave, MW, = Mach 
Number on cone surface, and @; = cone half-angle. Two widely 
These are 20° and 7.5°. 


In such a case, the author would 


‘ 


‘appreciably”’ 


different values of 6, are considered. 
An examination of the two tables indicates that “appreciably”’ 


is not very much. 


TABLE 1 


0, = 20° 
Us + Mo Ou M> Oy M, 
0.40 1.4 1.319 57.698 1.128 4.762% 0.976 
1.405 1.314 58.043 1.438 4.737 0.970 
0.50 1.4 1.661 44.206 ~ 1.429 6.737 1.291 
1.405 1.653 44.424 1.422 6.690 1.283 
0.60 1.4 2.141 35.981 1.797 9.251 1.677 
1.405 2.130 36.124 1.788 1” 9.196 1.667 
0.70 1.4 2.853 30.340 2.292 11.955 2.192 
1.405 2.839 30.435 2.280 11.903 2.178 
0.89 1.4 4.174 26.141 3.060 14.682 2.981 
1.405 4.154 26.199 3.042 14.640 2.963 
0.85 1.4 5.571 24.391 3.675 16.030 3.608 
1.405 5.546 24.433 3.653 15.997 3.586 
TABLE 2 
0, = 7.5° 
Us 1 Mo Ow M, Op M, 
0.50 1.4 1.366 47 .364 1.358 0.225 1.291 
1.405 1.358 47.728 1.350 0.219 1.283 
0.60 1.4 i Free 34.734 ter ogg 0.415 1.677 
1.405 1.761 34.962 1.747 0.403 1.667 
0.70 1.4 2.324 26 .082 2.292 0.806 2.192 
1.405 2.311 26.239 2.279 0.792 2.178 
0.80 1.4 3.198 19.329 3.109 1.623 2.981 
1.405 3.180 19.430 3.091 1.601 2.963 
0.90 1.4 5.133 13.506 4.776 3.340 4.617 
1.405 5.103 13. 560 4.749 3.316 4.588 
e 


Displacement and Rotation at the Junction of 
Two Cylinders of Unequal Thickness Due to 
Axisymmetrical Loadings 


G. D. Galletly 


Mechanical and Electrical Engineering Department, Shell 
Development Co., Emeryville, Calif. 


August 14, 1958 


i IS FREQUENTLY DESIRABLE to know the rotation and dis- 
placement at the junction of two cylinders of unequal thick- 
ness when axisymmetric loadings are applied around the junc- 


SPA es 
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tion. Such information is required, for example, if a transverse 
bulkhead is placed at the shell junction, and an analysis of bend- 
ing stresses is required. It is the purpose of the present note to 
obtain the M-@ relations, etc., for three axisymmetrical loadings of 
interest. These are 

(1) Uniformly distributed moments, M. 

(2) Uniformly distributed radial forces, F 

(8) Uniformly distributed moments and radial forces with 
radial displacement of the junction prevented 

The basic equations used are those relating the displacements 
and rotations at the end of a semi-infinite cylinder in terms of 
the edge moments and shears. These well-known equations 
may be written as 


w = (1/2Dh%) (AM — V (1) 
@ = (1/2Dd?) (2AM — V) (2) 
where 
w = edge radial displacement, + outwards 
6 = edge rotation, + counterclockwise 
M = edge bending moment 
V = edge shear force 
and D = Et®/12(1 — v?) (3) 
4 = 3(1 — 2p?) ‘R212 (4) 


where ¢ = thickness of shell and R = mean radius of shell 
Eqs. (1)-(4) are then applied to both the upper and lower 
shells, and the problems are solved by satisfying the equations of 


continuity and static equilibrium. The continuity relations for 
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Uniformly distributed radial forces acting around 
junction 


Fic. 2. 


all the problems are that the edge displacements and rotations 
of both cylinders be equal (for Case 3, w = 0). The equations 
of equilibrium for Cases 1 and 2 can be obtained from Figs. 1 and 2 
For Case 1, there are no external forces acting. Hence, the shear 
forces acting on the upper and lower cylinders must be equal and 
opposite; the external moment / must also be such that 


M = My, + M (5) 


where the subscripts u and / refer to the upper and lower cylinders, 
respectively. 

For Case 2, there are no external moments acting. 
the internal moments acting on the upper and lower cylinders 
must be equal and opposite: the relation between the internal 


Thus, 


shear forces V and the external force F is 
F= V, + V; (6) 


For Case 3, in which both external moments and forces exist, 
equations similar to both Eqs. (5) and (6) must be satisfied 

The analytical work is fairly straightforward, consisting 
mainly of manipulation of equations. In consequence, in the 
examples following, only final results are presented. The follow- 
ing notation will also be used for brevity: 


B = D.d.2/Diarv? = (R1/Ru): (tu/ti)? (7) 
y = du/Ar = (Riti/Ruatu)"? (8) 
¢ = (1 + B)? + 26[y + (1/y)] (9) 
Case (1). Uniformly Distributed Moments: 
6 = [(1 + By)/¢] (M/D~Ad (10) 
w = [(8 — 1)/t] (M/2D,,?) (11) 
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ne 


Case (2 Uniformly Distributed Radial Forces: 
6 = —[(8 — 1)/¢] (F/2DA2 12 
w= —)}(1 + (B/y)]/¢} (F/2DA3 13 
Case (3). Uniformly Distributed Moments and Radial Forces 

Radial Displacement at Junction Prevented: 

This problem can be solved by superposing Cases 1 and 2 
For no radial displacement at the junction, there follows from 
Eqs. (11) and (13): 

F/M = [(8 — 1)Ai]/[1 + (B/3 14) 

Inserting Eq. (14) into Eq. (2), and adding the result to Eq 
(10), one obtains 

6 = {1/[1 + (B/y)]} (M/2D~, 15 

Actually, Case 3 can be solved very simply without reference 
to Cases 1 and 2. 
there follows from Eq. (1) 


Vu = AML, 


This is because— for no radial displacement 


Vi = AM, 16) 


Insertion of the relations (16) into Eq (2), together with the 


use of Eq. (5), yields Eq. (15) immediately 


Errata—‘‘Free Molecular Flow Forces and Heat 
Transfer for an Infinite Circular Cylinder at 
Angle of Attack’’* 


L. Talbot 
Associate Professor, College of Engineering, University of 
California, Berkeley, Calif. 


October 14, 1958 


D* W. A. GustaFrson of Lockheed Aircraft Corp., Sunnyvale, 
Calif., has kindly called to my attention two errors which 
appeared in my article 

In Eq. (15) the correct expression for f(£) is that given, multi- 
plied by a factor +/z which was omitted. In Fig. 3 the labeling 
of curves f() and g() should be reversed 
* Readers’ Forum, Journal of the Aeronautical Sciences, Vol. 24, No. 6, 
pp. 458-459, June, 1957 


Pressure Fluctuations in a Blowdown Wind 
Tunnel 


James S. Murphy 


Design Specialist, Douglas Aircraft Co., 
Inc., El Segundo Division, El Segundo, Calif. 


August 9, 1958 


O™ OF THE PROBLEMS encountered in the operation of a blow- 
down wind tunnel is the occurrence of large-amplitude 
pressure fluctuations in the stilling chamber and test section. 
Preliminary measurements in the Douglas Trisonic 1-ft. wind 
tunnel indicated that the magnitude of these fluctuations was 
approximately 2 to 5 per cent of tunnel stagnation pressure for 
the test-section Mach Number range 0.2 < Mrs. < 18 A 
simple analysis of the flow in the stilling chamber (tunnel con- 
traction ratio at Mach 1 is 11.5) indicated that this fluctuation 
level was from 10 to 200 times the local dynamic pressure. 
Hence, it was concluded that the pressure fluctuations were not 
turbulence in the usual subsonic wind-tunnel sense—i.e., vor- 
ticity—since the upper limit for this type of fluctuation is of the 
order of the dynamic pressure. 

As is well known, two types of disturbances in addition to vor- 
ticity may exist in the flow of a compressible fluid—namely, 
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temperature or entropy spottiness and sound waves.! The air 
in the stilling chamber of the Douglas tunnel has the following 
history: 

(a) It is stored at a pressure of 125 psig and normal temper- 
ature (100°F.) in an 8,000-ft.* steel reservoir. 

(b) It is withdrawn from the reservoir during a blowdown 
period of approximately 40 sec., during which time the pressure 
and temperature of the reservoir air decreases while heat is added 
from the storage reservoir walls and a thermal mass at the down- 
stream end of the reservoir. 

(c) It passes over a 12-in. diameter butterfly valve whose posi- 
tion is controlled so that the air is throttled from storage pressure 
to a preselected, constant, stilling-chamber pressure. 

A certain amount of temperature spottiness is probably intro- 
duced by the various heat-transfer processes which occur up- 
stream of the stilling chamber. In addition, the flow in the 
neighborhood of the butterfly valve is sonic, and local supersonic 
regions exist with various oblique and normal shock waves which 
produce the desired throttling between the storage reservoir and 
stilling chamber. It is known that shock waves are excellent 
devices for generating noise as a result of the interaction of con- 
vected vorticity with the wave.” 

From microphone measurements of sound intensity outside 
the duct containing the butterfly valve, temperature-fluctuation 
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Muffler installation in stilling chamber. 


measurements in the stilling chamber using a hot-wire anemom- 
eter, and pressure fluctuations in the stilling chamber usiag a 
high-response, flush-diaphram type, strain-gage pressure trans- 
ducer, it was determined that the primary cause of the large pres- 
sure fluctuation level in the stilling chamber was a high-intensity 
sound field, superimposed on the mean flow and originating in 
the neighborhood of the control valve. Hence, it was felt that a 
good approach toward reduction of the fluctuation level would 
be to install sound- absorbent material in the form of a muffler in 
the stilling chamber. 

After the examination of several approaches to such a muffler, 
the design submitted by International Aerocoustics Corp. was 
chosen for installation in the tunnel between the entry duct 
and screen section, as shown in Fig. 1. This muffler consists of 
seven airfoil-shaped acoustic panels which span the 46-in. I.D. 
stilling chamber. 

Each panel is 4 in. wide and consists of a perforated galvanized- 
steel outer shell packed with bonded fiberglass. In order to pre- 
vent erosion of fiberglass through the holes in the shell, the 
fiberglass is enclosed in fiberglass-cloth bags which are covered 
with 100-mesh stainless-steel screening. The panels are on 8 in. 
centers. 
than in the entry duct because the lower velocity past the panels 
also reduces the possibility of fiberglass erosion. 

Measurements of the pressure-fluctuation level in the stilling 
chamber and test section before and after installation of the 
muffler are shown in Fig. 2. indicate that the 
muffler has reduced the fluctuation level by a factor of four or 
five, or by 12-14 decibels. The only penalty involved with this 
installation is a reduction of run time by about 1 sec. because of 


These data 


the increase of stilling chamber volume. 
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Transient Analysis of a Two-Chamber 
Blowdown Process With Critical Flow 


C. Forbes Dewey, Jr. 
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SYMBOLS 
m = mass, lb. 
y = ratio specific heats, dimensionless 
R = gas constant for flow medium, (ft. lb.)/(Ib. °R.) 
T = temperature,°R. 


The muffler was installed in the stilling chamber rather } 
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eo NOTE PRESENTS a simplified analysis of the transient con- 
ditions existing in a supply tank of given initial pressure and 
a second tank or receiver which in turn discharges to the atmos- 
phere. Although the analysis excludes the effects of complex 
expansion shocks and heat transfer, the results are readily appli- 
cable as a first approximation. Results for the limiting case of 
small transient times are derived in closed form 

Fig. 1 shows the system under consideration. The supply tank 
is connected to the second chamber by a nozzle of area a. 
Chamber 2 
tem through a throat of area a». 


discharges to the atmosphere or a low-pressure sys- 
The process is assumed to be 
adiabatic and the flow through areas a; and ad is considered to be 
given by the isentropic steady-state relation for critical flow 


dm/dt = —V vg RT pa[l + (y — 1)/2] 7%" I/2y-1) (1) 
where 7 and p are the instantaneous stagnation values of tem- 
perature and pressure calculated from the energy equation for a 
perfect adiabatic gas. This implies that the time required to 
reach a receiver pressure giving critical flow through area a2 is 
negligible compared to the total transient time, and that the ratio 
pi/p2) is large enough to maintain critical flow through a, dur- 
ing the entire process. 

For the first tank, the well-known adiabatic relations between 


state variables apply. They may be written in the form 


T:/Tn = {1 + (vy — 1)/2] VTu/Bit}-? = 


(7,/m)% = ¥ {3} 


(pi/Po \(¥ —1) 
/a;*yRg, «+ =1,2 (8) 


Vi? [1 + (y — 1I)/Q]OFYD/O-) 


Considering the second chamber, the continuity and energy 


equations are written 


dm. = —dm, — mov T>2/B2 dt (4) 


(1/y)Todmz + (1/y)medT2 = —Tidm, — T2V T2/B2 modt (5) 


Combination of the above equations gives the mass and tem- 
perature in the second chamber: 
dT, = }(yTumo/m2)V Tu/B: (1 + (vy — 1) X 


— (37 —-1)/(y-1) yTon Ti. Bo} dt — 


T2dm2/me 


VT. Bi t 2] 
(6) 
(y — 1) VWTa/fr t/2) 7 VF Dy - DS 


dmz = jmuV To /Bi [1 + 


meV To Bo} dt (7) 

Eqs. (6) and (7) now represent two different equations con- 
taining the three variables me, 72, and ¢. Under special condi- 
tions, these may be simplified and solved directly; in any case, 
the evaluation may be done numerically. 

If the time required to establish the desired operating pressure 
in tank 2 is small, or if the supply tank is very large, Eqs. (6) and 
(7) may be reduced by assuming that the pressure and temper- 
ature in the supply tank are constant. For the case of equal 
initial temperatures, the equations may be nondimensionalized 


with respect to the initial values to give 


dT,* = (1/m2*)V Ti2/B> ((po/pi2) (ai/a2) (y — T2*) — 
(y — 1)T2*mo*V T2*] dt (8) 
dmz* = V T.2/B2 ((pu/Po2) (a1/a2) — ms*/ T2*)dt (9) 


During the first part of the transient period, the temperature 
y; as m.* increases, (d7.*/dt) changes sign 
If (poi/poz) (a:/az) > m2*V/ T2*, 


will approach 72* = 
T2* falls off toward 1. 
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Fic. 1. Time required to reach chamber pressure 


representing 7)* by a constant will yield a reasonable approxima 


tion to the continuity equation and the temperature 7)* may be 
adequately determined from the relation 
T-* = (7 + 62/2) - bN 7 + b2 1| 10 
(10) 
b= p2*( y — 1)/(poi/ poz) (ai/a2) f 
Of particular interest is the time required to reach a given 


pressure in the second tank; this is obtained by integration and is 
written 


(1/V T2* V Te B2) log 1 (por Por) (i /d2) — 
V T:*] [(por/Po2) (ai/a2) — po*/V T2*}} 


i= 
(11) 
where 7>* is given by Eq. (10). 

Substituting the values shown in the figure into Eqs. (10) and 
(11), it is apparent that the time required to reach a given oper- 
If the 
throat area ad» is fixed by design considerations, the curve will 


ating pressure is quite sensitive to the area ratio (a@;/a 


immediately show what supply area a, is necessary to establish 
the required pressure within a given time period 

After preliminary estimates have been made of the system 
design dimensions using calculations based on the assumption of 
small time, Eqs. (6) and (7) may be easily solved numerically to 
As a final 


refinement, the effects of heat transfer may be accounted for as 


obtain a more adequate description of the operation. 
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Theoretical Pressure Distribution on a 
Hemisphere-Cylinder Combination; 
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i RECENT YEARS great use has been made of approximate 
methods for the determination of the pressure distribution 
on blunt-nosed bodies and afterbodies at high Mach Numbers 
For quasi-spherical bodies it has been suggested that modified 
Newtonian theory in combination with a Prandtl-Meyer expan- 

+ This investigation was under Contract AF 49(638)-217, 
sponsored by the USAF through the Air Force Office of Scientific Research 
Air Research and Development Command, in connection with reference 4 

* The author wishes to express his gratitude to Dr. Antonio Ferri and 
Dr. Roberto Vaglio-Laurin for their helpful suggestions 
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sion be used on the nose portion, the two laws being matched 50; 2 a GE | 

at the point where the pressure gradients are equal. No simple 45: tot MODIFIED NEWTONIAN| 
approximation, however, has been found for flat-nosed bodies. 

As for the pressure distribution on the afterbody, the blast-wave 40; 
analogy*® has been suggested for general nose shapes but par- 35+ 
ticular afterbody profiles. 

The purpose of the present note is to compare these approxi- 30) ’ CHARACTERISTICS | 
mate estimates with a more accurate determination of the flow 25. } 
field about a hemisphere-cylinder in an ideal gas flow at M = 20. . \h J 7 PRANDTL-MEYER 
It was felt that since experimental investigations in air at this 20) rN 74 
Mach Number are scarce and very difficult to obtain, the com- IS} 


4 rie j BLAST WAVE ORIGIN ONE RADIUS UP- 
parison would be of interest. The basis of comparison is the 10! 


ID APPROXIMATION 
flow field as it results from a numerical integration of the exact See ==BUAST lig SECOND 
equations governing the motion of the ideal fluid. The integra- 5) rT __ APPROXIMATION) 
tion in the elliptic region has been performed by the method of o.— = SS ae oD Neos ee J 
reference 1 and has been continued in the hyperbolic region, 0 20 40 60 na 0 20 40 160 180 
outside the domain of influence of the sonic line, by the standard R, 
characteristic method? for rotational axisymmetric flow. The ‘ ’ oe . ; , 

ees : 3. 2. Pressure distribution on a hemisphere-cylinder af 
results are presented in Figs. 1, 2, and 3. It appears that the Me = 20 (cylindrical part). , 
modified Newtonian-plus-Prandtl-Meyer flow values are in rea- 
sonably good agreement with the more accurate results on the 
nose portion of the body. As for the cylindrical afterbody, it 
appears that the blast-wave analogy, taken to a second approxi- 
mation, yields the proper trend for the shock-wave shape and 
for the pressure decay; however, the predicted pressures are siz- 
ably higher than those determined by the characteristic calcula- 
tion. Better agreement could be obtained by choosing a suit- 
able origin for the blast. This shift in coordinate axis would be 
valid since the solution in question is an asymptotic one. The 
present data would suggest an effective origin located approxi- t 
mately five radii upstream of the nose; such a shift, however, is , | es 
much larger than that suggested by Lees* on the basis of a com- 
parison with experimental data at JJ = 7.7. Therefore, it ap- 
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pears that the proper choice of the origin depends on the Mach 
Number considered, and a general rule is not yet available. 
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GIANNINI PLASMADYNE CORPORATION 

THE B. F. GOODRICH COMPANY 

GOODYEAR AIRCRAFT CORPORATION 

GRUMMAN AIRCRAFT ENGINEERING CORPORATION 

HYDRO-AIRE, INC. 

INSURANCE COMPANY OF NORTH AMERICA COM- 
PANIES 


INTERNATIONAL BUSINESS MACHINES CORPORATION 
THE INTERNATIONAL NICKEL COMPANY, INC. 


ITT LABORATORIES, DIVISION OF INTERNATIONAL 
TELEPHONE AND TELEGRAPH CORPORATION 


JANITROL AIRCRAFT DIVISION, SURFACE COMBUS- 
TION CORPORATION 


JOHNS-MANVILLE SALES CORPORATION 
EARLE M. JORGENSEN CO. 
KELSEY-HAYES COMPANY 

WALTER KIDDE & COMPANY. INC, 
KOLLSMAN INSTRUMENT CORPORATION 


KOPPERS COMPANY, INC. 
SOUND CONTROL DEPARTMENT 


LAVELLE AIRCRAFT CORPORATION 
LEAR, INCORPORATED 

Cc. W. LEMMERMAN, INC. 
LIBRASCOPE, INCORPORATED 

THE LIQUIDOMETER CORPORATION 
LOCKHEED AIRCRAFT CORPORATION 


LOEWY-HYDROPRESS DIVISION OF BALDWIN-LIMA- 
HAMILTON CORPORATION 


LONGINES-WITTNAUER WATCH COMPANY, INC. 
MARQUARDT AIRCRAFT COMPANY 

THE MARTIN COMPANY 

MARTIN STEEL PRODUCTS CORPORATION 
McDONNELL AIRCRAFT CORPORATION 
MELETRON CORPORATION 

MID-CONTINENT MANUFACTURING, INC, 
MINNEAPOLIS-HONEYWELL REGULATOR COMPANY 
NATIONAL CREDIT OFFICE, INC. 

NORTH AMERICAN AVIATION, INC. 

NORTHROP AIRCRAFT, INC. 


NUCLEAR DEVELOPMENT CORPORATION OF A 
PAN AMERICAN WORLD AIRWAYS, INC. 
THE RALPH M, PARSONS COMPANY 


PESCO PRODUCTS DIVISION, BORG-WA 
PORATION 


PHILLIPS PETROLEUM COMPANY 

PIASECKI AIRCRAFT CORPORATION 

THE RAMO-WOOLDRIDGE CORPORATION 
REPUBLIC AVIATION CORPORATION 

ROHR AIRCRAFT CORPORATION 

PAUL ROSENBERG ASSOCIATES 

RYAN AERONAUTICAL COMPANY 
SANDBERG-SERRELL CORPORATION 
SCHRILLO AERO TOOL ENGINEERING COMPAR 
SHAFER BEARING DIVISION, CHAIN BELT CO! 
SHELL OIL COMPANY 

SIMMONDS AEROCESSORIES, INC, 

SOCONY MOBIL OIL COMPANY, INC. 

SOLAR AIRCRAFT COMPANY 

SOUTHWEST PRODUCTS CO, 

R. DIXON SPEAS ASSOCIATES 


SPERRY GYROSCOPE COMPANY DIVISION OF 
RAND CORPORATION 


STANDARD OIL COMPANY OF CALIFORNTA 
STANDARD OIL COMPANY (INDIANA) 


STANDARD-THOMSON CORPORATION 
CLIFFORD MANUFACTURING COMPANY DI 


STANLEY AVIATION CORPORATION 
STROUKOFF AIRCRAFT CORPORATION 
THIOKOL CHEMICAL CORPORATION 
THOMPSON RAMO WOOLDRIDGE INC, 
TOOLKO ENGINEERING COMPANY 

TRANS WORLD AIRLINES, INC. 

TURBO PRODUCTS, INC. 

UNION CARBIDE CORPORATION 

UNITED AIR LINES, INC. 

UNITED AIRCRAFT CORPORATION 

UNITED STATES AVIATION UNDERWRITERS, IN 
VARD, INC. 

VERTOL AIRCRAFT CORPORATION 

VICKERS INCORPORATED 

VITRO CORPORATION OF AMERICA 

WESTERN GEAR CORPORATION 
WESTINGHOUSE ELECTRIC CORPORATION 
— INSTRUMENTS DIVISION OF DAYST 


WYMAN-GORDON COMPANY 
YOUNG RADIATOR COMPANY 























